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FOREWORD 


This  document  presents  the  models,  analysis 
and  results  obtained  for  the  National  Weather  Service 
under  Contract  No.  NA79SAC00668 ,  "Applications  of 
Kalman  Filtering  and  Maximum  Likelihood  Parameter 
Identification  to  Hydrologic  Forecasting". 

A  computer  program,  "SUBROUTINE  REDO-UHG" , 
and  a  report  entitled  "Reduced  Order  Unit  Hydrograph 
Program  Documentation"  providing  information  on  the 
design  and  use  of  the  program  have  been  previously 
delivered  to  the  National  Weather  Service  as  part  of 
the  same  contract. 

This  study  has  benefited  from  several 
conversations  with  E.A.  Anderson,  E.R.  Johnson  and 
G.  Smith  of  the  Hydrologic  Research  Laboratory  of 
the  National  Weather  Service  and  with  Professor  R.L. 
Bras  of  MIT  and  his  assistants.  Their  cooperation 
is  gratefully  acknowledged. 
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ABSTRACT 


The  applications  of  the  canonical  variate, 
Kalman  filtering  and  maximum  likelihood  parameter 
identification  techniques  to  the  requirements  of  the 
National  Weather  Service  in  river  flow  forecasting 
are  investigated. 

State  space  reduced-order  models  for  unit 
hydrographs  are  obtained  with  the  use  of  canonical 
variate  methods.  A  complete  state-space  model  for  a 
catchment  consisting  of  the  Sacramento  model  as  the 
soil  moisture  system  and  the  basin's  unit  hydrograph 
as  the  channel  routing  system  is  constructed.  This 
model  is  used  in  the  design  of  extended  Kalman  fil¬ 
ters  for  the  prediction  of  the  channel  discharge  and 
the  state  of  the  system,  and  also  in  the  design  of 
an  algorithm  for  the  identification  of  catchment 
model  parameters  through  the  use  of  maximum  likeli¬ 
hood  techniques.  The  performance  of  the  algorithms 
is  demonstrated  with  synthetic  data  generated  with 
the  models  for  the  Bird  Creek  and  White  River  basins 
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1 .  INTRODUCTION 

1 .  1  BACKGROUND 

The  National  Weather  Service  (NWS)  has  the  responsi¬ 
bility  for  hydrologic  forecasting  in  the  United  States.  This 
responsibility  includes  the  production  of  both  flood  warnings 
and  stream-flow  forecasts.  Accurate  and  timely  flood  warnings 
are  required  for  a  wide  variety  of  flood  classes  including 
flash  floods  as  well  as  floods  of  longer  duration.  Stream-flow 
forecasts  are  required  for  diverse  applications  including  the 
planning  of  irrigation,  the  prediction  of  available  hydroelec¬ 
tric  power,  the  maintenance  of  water  quality  standards,  and 
the  planning  of  river  navigation. 

There  is  a  continuing  need  for  new  techniques  useful 
for  creating  more  accurate  and  cost-effective  flood  warning 
and  stream-flcw  predictions.  It  is  highly  desirable  to  in¬ 
crease  the  amount  of  automation  used  in  the  creation  of  hydro- 
logic  forecasts  and  to  be  able  to  take  advantage  of  newly- 
developing  advances  in  computer  and  communications  capabilities 
and  in  computational  and  algorithmic  techniques. 

1.2  STUDY  OBJECTIVES 

The  objective  of  this  study  was  to  investigate  the 
application  of  Kalman  filtering,  canonical  variate  and  maximum 
likelihood  parameter  identification  techniques  to  the  require¬ 
ments  of  the  National  Weather  Service  in  improving  hydrologic 
forecasting.  The  work  was  organized  into  three  principal  tasks 
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Task  1  ISSUES  IN  FILTER  DESIGN 

Under  this  task,  the  differential  equations  corre¬ 
sponding  to  the  Sacramento  Soil  Moisture  model 
were  derived  and  several  theoretical  modifications 
were  developed.  The  equations  of  the  soil  moisture 
model  were  combined  with  a  reduced-order  state- 
space  model  for  a  unit  hydrograph  and  an  extended 
Kalman  filter  that  produces  six-hour  lead  forecasts 
of  a  basin's  discharge  was  designed  and  implemented. 
The  details  of  the  analysis  and  a  set  of  results 
obtained  with  the  SSM  model  parameters  of  the 
Bird  Creek  and  White  River  drainage  basins  are 
included  in  this  report. 

Task  2  STATE-SPACE  MODEL  DEVELOPMENT  FOR  UNIT  HYDROGRAPHS 

Under  this  task,  a  computer  program  has  been 
developed  for  applying  the  canonical  variate 
technique  to  the  development  of  di sc  re t e- t ime 
reduced-order  state-space  models  for  the  approx¬ 
imation  of  unit  hydrographs.  The  computer  pro¬ 
gram  (Subroutine  REDO-UHG)  accompanied  with  sup¬ 
porting  documentation  (Ref.  1)  has  been  delivered 
to  NOAA/NWS  for  use  as  an  operation  in  the  Version 
5.0  NWSRFS  Forecast  Component.  The  principles 
on  which  the  design  of  the  computer  program  was 
based  are  described  in  this  report  together  with 
some  examples  of  their  application. 

Task  3  PARAMETER  IDENTIFICATION  FOR  CATCHMENT  MODELING 

Under  this  task,  an  initial  investigation  of  the 
applications  of  the  technique  of  maximum  likeli¬ 
hood  parameter  identification  to  the  problem  of 
catchment  calibration  has  been  performed.  Pa¬ 
rameter  estimation  algorithms  appropriate  to  the 
catchment  model  of  the  National  Weather  Service 
have  been  developed  and  tested  with  simulated 
data  to  determine  parameter  estimation  error, 
parameter  i den t i f i abi 1 i t y  and  numerical  behavior 
of  the  algorithms. 


Figure  1.2-1  presents  the  hierarchical  structure  of 
the  outputs  of  the  three  tasks  described  above.  The  state- 
space  models  of  un i t - hyd rographs  obtained  with  the  canonical 
variate  technique  under  Task  2  have  direct  application  to  t  h  t 
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simplification  of  the  processing  of  channel - inf low  time  series. 
These  models  for  unit  hydrographs  are  used  in  Task  1  to  create 
complete  state-space  models  for  catchments  which,  in  turn, 
form  the  basis  for  the  design  and  operation  of  Kalman  filters 
for  the  prediction  of  a  basin's  channel  discharge  and  for  the 
estimation  of  the  state  of  the  system. 

The  maximum  likelihood  parameter  identification  proce¬ 
dure  of  Task  3  is  an  iterative  algorithm  which  uses  a  Kalman 
filter  as  one  of  its  main  constituents.  The  filter  innovations 
(differences  between  actual  and  predicted  measurements)  are 
used  in  the  evaluation  of  the  likelihood  function.  An  optimi¬ 
zation  procedure  utilizes  the  values  of  the  likelihood  function 
and  its  functional  gradient  computed  by  symbolically  differen¬ 
tiating  the  operations  of  the  filter  to  evaluate  a  vector  param¬ 
eter  increment  in  the  direction  that  maximizes  the  likelihood 
function.  The  catchment  model  parameters  are  modified  and  the 
process  repeated  until  a  convergence  criterion  is  satisfied. 


1.3  REPORT  ORGANIZATION 

The  organization  of  this  report  is  as  follows:  Sec¬ 
tion  2  deals  with  the  catchment  model  and  the  design  of  ex¬ 
tended  Kalman  filters  for  the  prediction  of  a  basin's  channel 
discharge.  Section  3  presents  the  canonical  variate  technique 
and  its  application  to  the  synthesis  of  reduced-order  state- 
space  models  for  unit  hydrographs.  Section  A  describes  the 
algorithms  for  maximum  likelihood  identification  of  the  param¬ 
eters  of  the  catchment  model  and  presents  some  examples  of 
their  application.  Section  5  contains  a  review  of  the  infor¬ 
mation  contained  in  this  report.  An  appendix  containing  some- 
technical  considerations  on  the  constraint  on  the  ratio  of 
free  to  tension  water  for  the  upper  zone  of  the  soil  moisture 
model  is  also  included. 
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For  the  sake  of  completeness,  some  of  the  material 
contained  in  previous  progress  reports  (Refs.  2  through  5)  is 
also  included  in  this  report. 
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2. 


ISSUES  IN  FILTER  DESIGN 


This  section  presents  the  details  of  the  analysis  and 
synthesis  of  an  Extended  Kalman  Filter  (EKF)  for  improving  the 
estimation  of  the  state  and  channel  discharge  of  a  basin.  The 
basin's  model  considered  in  this  study  consists  of  the  Sacramento 
model  as  the  soil  moisture  system  and  a  unit  hydrograph  as  the 
channel  routing  system. 

Pioneering  work  on  the  applications  of  Kalman  filter¬ 
ing  to  the  National  Weather  Service's  river  flow  forecasting 
system  is  described  in  Ref.  6.  The  present  work  differs  con¬ 
ceptually  from  that  of  Ref.  6  in  that  the  soil  moisture  proc¬ 
ess  is  viewed  as  occurring  in  conti nuous- t ime  rather  than  in 
di sc rete- t ime  as  in  Ref.  6."  In  addition,  channel  routing  is 
modeled  by  the  basin's  unit  hydrograph  while  in  Ref.  6  a  linear 
reservoir  with  variable  outflow  rate  is  used. 

The  Kalman  filtering  formulation  requires  that  the 
system  be  modeled  in  state-space  form.  The  state-space  de¬ 
ferential  equations  of  the  Sacramento  Soil  Moisture  (SS.M)  model 
were  derived  in  Refs.  2,  3,  and  A.  To  complement  these  equa¬ 
tions,  channel  routing  was  modeled  by  a  suitable  order  state- 
space  model  approximation  to  the  basin's  unit  hydrograph  ob¬ 
tained  with  the  methods  described  in  Chapter  3  of  this  report. 

The  basic  idea  behind  the  EKF  formalism  is  to  approxi¬ 
mate  the  system  behavior,  for  a  short  time  interval,  by  the 


*A  comparison  of  the  state-equations  of  Ref.  (>  and  those  used 
in  the  present  study  is  given  in  Ref.  2.  Other  differences 
between  the  two  studies  are  also  noted  in  Ref.  2. 
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linearized  form  of  the  state-equations  about  the  operating 
point  at  the  beginning  of  the  interval.  The  optimal  filter 
for  the  resulting  linear  system  is  used  to  propagate  the  state 
estimates  and  their  covariance  matrix  during  the  time  interval 
under  consideration.  A  detailed  treatment  of  extended  Kalman 
filtering  is  given  in  Ref.  7. 

In  order  to  test  the  performance  of  the  filter  and 
the  consistency  of  the  formulation,  synthetic  data  was  gener¬ 
ated  using  the  model  in  the  simulation  mode.  These  data  were 
then  used  as  input  to  the  filter,  and  the  state  and  discharge 
estimates  so  obtained  were  compared  to  the  truth  values  previ¬ 
ously  generated.  The  SSM  model  parameters  and  unit  hydrographs 
used  in  these  tests  were  those  of  the  National  Weather  Service 
River  Forecasting  System  (NWSRFS)  calibration  of  the  Bird  Creek 
and  White  River  drainage  basins. 

This  chapter  is  organized  as  follows:  Section  2.1 
discusses  the  overall  model  structure  and  the  interface  be¬ 
tween  the  soil  moisture  model  and  the  system  associated  with 
the  unit  hydrograph.  Section  2.2  presents  the  state-equations 
of  the  cont inuous- t ime  part  of  the  model  including  the  SSM 
model  and  a  simple  precipitation  model.  Section  2.3  discusses 
the  discrete-t ime  equations  associated  with  the  unit  hydrograph 
system.  Sections  2. A  and  2.5  describe  in  detail  the  operation 
of  the  model  in  the  simulation  and  filtering  modes,  respectively. 
Section  2.6  presents  a  collection  of  representative  results 
obtained  with  the  techniques  described  in  previous  subsections. 

2.1  MODEL  STRUCTURE 

Traditionally,  NWS  has  used  deterministic  models  in 
forecasting  river  flows  based  on  meteorological  data.  Thus, 
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necessarily,  the  intrinsic  rate  at  which  the  soil  moisture 
model  operates  has  been  the  same  as  the  rate  at  which  precipi¬ 
tation  data  are  collected.  When  a  stochastic  conceptual  model 
is  used,  however,  the  input  rate  to  the  soil  moisture  model 
need  not  be  the  same  as  the  rate  at  which  accumulated  channel 
inflow  estimates  are  produced.  In  fact,  with  a  con t i nuous - t ime 
stochastic  soil  moisture  model,  the  best  estimate  of  the  channel  - 
inflow  rate  given  all  previous  measurements  of  precipitation 
and  channel  discharge  can  be  computed  at  any  time. 

Thus,  even  though  the  SSM  model  is  a  discrete- t ime 
model,  it  was  necessary  to  derive  a  cont i nuous - t i me  model 
whose  discretized  version  was  in  congruence  with  the  SSM  model. 

In  addition  to  the  advantage  of  being  able  to  compute  state- 
estimates  at  any  time,  there  were  two  other  reasons  for  model¬ 
ing  the  soil  moisture  process  in  cont inuous- t ime .  First,  the 
physical  processes,  in  themselves,  take  place  continuously  in 
time.  For  example,  the  effects  of  a  severe  storm  of  short 
duration  cannot  be  properly  modeled  using  a  predetermined 
equally  spaced  sequence  of  times.  Evidence  for  the  need  of 
considering  the  dynamic  behavior  of  the  system  in  continuous¬ 
time  is  found  in  the  LAND  subroutine.  There,  depending  on  the 
availability  of  free  water  in  the  upper  zone,  the  basic  time 
interval  is  partitioned  into  a  number  of  subintervals  for  the 
computation  of  the  percolation  function.  This  computation 
determines  the  distribution  of  water  to  the  lower  zone,  the 
amount  of  surface  runoff,  etc.  Secondly,  the  threshold  values 
associated  with  many  of  the  variables  can  be  attained  at  times 
which,  in  general,  do  not  coincide  with  the  endpoints  of  arbi¬ 
trarily  chosen  time  intervals.  These  thresholds  determine 
when  the  system  switches  from  one  mode  of  operation  to  another 
and  are  of  fundamental  importance  in  the  analysis. 
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Channel  routing,  on  the  other  hand,  is  modeled  as  a 
di sc rete- t ime  system.  A  k-hour  unit  hydrograph  yields  instan¬ 
taneous  discharge  rates  when  given  an  input  sequence  consisting 
of  k-hourly  accumulated  channel  inflow  values.  For  most  of 
the  unit  hydrographs  used  by  NWS  the  input  sample  rate  is  equal 
to  the  output  sample  rate;  i.e.,  estimates  of  the  discharge 
rate  are  produced  at  intevals  of  k  hours.  All  hydrographs 
considered  in  this  chapter  belong  to  this  category. 

In  the  present  analysis  measurements  of  accumulated 
precipitation  are  assumed  to  occur  at  an  £k-hourly  rate  and 
observations  of  mean  discharge  occur  every  m£k  hours.  For 
example,  the  situation  in  which  k=6,  fi=l,  m=l  represents  the 
case  where  continuous  estimation  of  the  state  of  the  system  is 
performed  given  6-hour  measurements  of  accumulated  precipita¬ 
tion  and  instantaneous  discharge;  k=6,  fi=l,  m=4  corresponds  to 
6-hour  measurements  of  precipitation  and  daily  observations  of 
mean  discharge.  Figure  2.1-1  summarizes  the  above  convention. 

In  order  to  combine  the  state-space  model  of  the  soil 
moisture  accounting  procedure,  which  yields  cont i nuous- t ime 
estimates  of  the  channel  inflow  rate,  with  the  unit  hydrograph 
system,  which  requires  k-hourly  accumulated  channel  inflow  at 
its  input,  it  is  necessary  to  introduce  an  additional  state. 

The  role  of  this  state  is  to  integrate  the  channel  inflow  rate 
for  periods  of  k  hours.  A  schematic  diagram  depicting  the 
interrelation  between  the  different  components  of  the  model  is 
given  in  Fig.  2.1-2.  The  additional  state  mentioned  above  can 
be  visualized  as  a  reservoir,  labeled  Channel  Inflow  Accumu¬ 
lator  in  Fig.  2.1-2,  whose  contents  are  dumped  into  the  unit 
hydrograph  system  every  k  hours. 


‘'■In  some  instances,  NWS  uses  unit  hvdrographs  for  which  the 
output  rate  is  higher  than  the  input  rate  (see  Section  1 . U  .  I )  . 
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Figure  2.1-1  Convention  on  Measurements  and  Unit 
Hydrograph  Operation  Rate 
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Figure  2.1-2  Basin  Model  Components 
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Two  more  accumulator  states  are  needed  to  account  for 
the  memory-type  measurements.  The  rainfall  accumulator  in  the 
upper  right  corner  of  Fig.  2.1-2  integrates  the  rainfall  rate 
for  periods  of  £k  hours  and  the  discharge  accumulator  in  the 
lower  left  corner  of  Fig.  2.1-2  adds  up  the  m£  discharge  rates 
whose  average  yields  the  m£k-hour  mean  discharge. 

The  basin's  state-space  model  contains  10+M  states 
where  M  is  the  order  of  the  state-space  model  approximation  to 
the  unit  hydrograph  system.  A  description  of  the  states  is 
given  in  Table  2.1-1.  It  was  shown  in  Ref.  2  that  six  states 
suffice  to  represent  the  SSM  model.  The  first  six  states  in 
Table  2.1-1  correspond  to  the  SSM  model.  Their  equivalents  in 
the  LAND  subroutine  are  indicated  in  parentheses  in  the  table. 
States  7,  9,  and  10+M  are  associated  with  the  accumulators 
previously  described.  State  8  provides  the  basis  for  the  rain¬ 
fall  model  which  is  presented  in  detail  in  Section  2.2.2.  The 
remaining  states  (10  through  9+M)  correspond  to  the  unit  hydro¬ 
graph  system. 

It  is  convenient  to  partition  the  state-vector,  x.  as 

‘■(t) 

where  xc  and  x^  stand  for  the  first  9  and  last  M+l  components 
of  x,  respectively.  The  subvector  xc  evolves  continuously  in 
time  while  x^  changes  only  at  times  which  are  multiples  of  the 
hydrograph  rate,  k.  These  times  (vk;  v=0,l,...)  are  referred 
to  as  critical  times  in  the  sequel  . 


The  operation  of  the  filter  can  be  described  in  general 
terms  as  follows.  Between  critical  times  there  are  n<>  measure- 
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TABLE  2.1-1 

BASIN  MODEL  STATE  DESCRIPTION 


STATE 

DESCRIPTION 

1 

Upper-zone  tension-water  content  (UZTWC) 

2 

Upper-zone  free-water  content  (UZFWC) 

3 

Lower-zone  tens  ion -wa t er  content  ( LZTWC ) 

A 

Lower-zone  primary  free-water  content  (LZFPC) 

5 

Lower-zone  supplementary  free-water  content 
( LZFSC ) 

6 

Excess  of  the  additional  impervious  storage 
over  the  upper-zone  tension-water  content 
( AD 1MC- UZTWC) 

7 

Channe 1 - in f low  accumulator  content 

8 

Rainfall  generator  model 

9 

Rainfall  accumulator  content 

10 

. 

/Unit  hvdrograph  model 

9+M 

) 

1 0+M 

Chanrel  discharge  accumulator 

merits.  The  continuous  part  of  the  state  estimate,  x  ,  is 
propagated  in  accordance  with  the  associated  differential  equa¬ 
tions,  but  the  discrete  part,  x^  ,  remains  unchanged.  When  a 
critical  time  is  reached,  the  channel  routing  portion  of  the 
model  is  updated  taking  into  account  the  accumulated  channel- 
inflow  during  the  past  k  hours.  At  this  point,  if  there  are 
any  measurements,  the  Kalman  gains  are  computed  and  used  to 
modify  all  the  state  values  by  incorporating  optimally  all  of 


x  is  the  estimate  of  x. 
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the  information  contained  in  the  measurements.  Finally,  the 
contents  of  the  appropriate  accumulators  are  reset  to  zero  and 
propagation  of  the  continuous  part  of  the  state  for  the  next  k 
hours  begins. 

It  is  convenient  to  introduce  a  notation  of  super¬ 
scripts  to  be  affixed  to  the  critical  times  in  order  to  dis¬ 
tinguish  the  various  state  values  computed  at  these  times.  If 
t  is  a  critical  time,  x(t  )  stands  for  the  state  vector  immedi¬ 
ately  after  propagation  of  the  continuous  part  for  the  last  k 
hours  has  been  completed,  x(t°)  represents  the  state  after  the 
discrete  transition  for  the  unit  hydrograph  part  of  the  model 
has  taken  place,  x(t+)  denotes  the  value  of  the  state  following 
a  Kalman  update  and  x(t  )  is  the  state  after  resetting  any 
accumulator  to  zero.  Table  2.1-2  summarizes  the  behavior  of 
the  state  values  for  a  full  cycle  of  operation  of  the  model 
(m£k  hours)  using  symbolically  the  notation  introduced  above. 
For  example,  the  notation  -  =  0  /  +  ^  r  means  x(t  )  =  x(t°) 

/  x(t+)  /  x(tr) . 

Thus,  the  transition  from  t  to  t°  only  affects  the 
discrete  states  of  the  model  (last  two  columns  in  Table  2.1-2). 
The  transition  from  t°  to  t+  corresponds  to  a  Kalman  update. 
Since  all  state  estimates  can  be  expected  to  improve  following 
a  Kalman  update.  Table  2.1-2  indicates  0  /  +  for  all  states  at 
the  times  at  which  there  is  at  least  one  measurement .  The 

•f 

transition  from  t  to  t  only  affects  the  accumulator  states. 
State  seven,  the  channel-inflow  accumulator,  is  reset  to  zero 
at  all  critical  times.  State  nine,  the  rainfall  accumulator, 
and  state  10+M,  the  discharge  accumulator,  are  set  to  zero 
following  a  rainfall  or  discharge  measurement,  respectively. 

The  arrows  in  the  last  two  columns  of  Table  2.1-2  indicate 
that  the  discrete  part  of  the  model  does  not  vary  between 
critical  times,  i.e.,  if  t.  and  t..,  are  two  successive  criti- 

l  j  +1 

cal  times,  then 
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W  =  ^d(ti+l) 


(2.1-2) 


2.2  MODEL  EQUATIONS  -  CONTINUOUS  TIME  COMPONENTS 

The  state-equations  for  the  complete  basin  model  are 
of  the  form  (v=0,l,...) 

xc(t)  =  Fc(xc,t)  +  GcC(t)  (2.2-1) 

5dl(vk)°J  =  A  xdl(vk)~)  +  B  xcl(vk)"J  (2.2-2) 

where  £c(x  , t)  is  a  nonlinear  t ime-varying  vector  function  of 
the  continuous  part  of  the  state,  the  time  dependency  being 
through  the  potential  evapot  ranspi  ra  t  i  on  demand,  u^(t).  C>c  is 
a  9><1  matrix  and  £  is  a  scalar  gaussian  white  noise  input  that 
drives  the  rainfall  model  (see  Section  2-2.2).  Thus,  the  only 
nonzero  entry  in  Gc  appears  in  row  eight.  The  factors  A  and  B 
are  (M+l)x(M+l)  and  (M+l)x9  constant  matrices,  respectively. 

In  addition  to  Eqs.  2.2-1  and  2.2-2,  the  accumula¬ 
tors'  contents  are  reset  to  zero  at  the  appropriate  times 

x7( (n)r|  =  0  ;  n  =  0  ,  k  ,  2k _  (2.2-2) 

x9((n)rl  =  0  ;  n  =  0,£k,2£k,  .  .  .  (2.2-9) 

X 1 0+M  I  (  n  )  r  I  =  0  ;  n=0  ,m£  k  ,  2m£  k  ,  .  .  .  (2.2-r)) 

Changes  in  the  discrete  part  of  the  state  vector,  x^, 
occur  only  at  the  critical  times.  If  vk  <  t  <  (v  +  1  )k  the  state- 
equations  for  xd  can  be  thought  of  to  be 

xd(t)  =  0  (2.2-6) 
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Let  and  y 2  be  the  measured  values  of  precipitation 
and  mean  discharge,  respectively.  The  measurement  equations  are 

y  j ( v  £  k )  =  x9!(vfik)°]  +  ^  2 ( v£  )  (2.2-7) 

y 2 ( v m£ k )  =  ^  x : Q+M [ ( vm£ k ) ° 1  +  42(vm£)  (2.2-8) 

where  and  ^2  rePresent  the  errors  in  the  measurements  mod¬ 
eled  as  white  noise  sequences. 

This  section  presents  the  individual  differential 
equations  associated  with  Eq .  2.2-1.  It  is  organized  as  fol¬ 
lows:  Subsection  2.2.1  treats  the  SSM  model  equations.  Sub¬ 

section  2.2.2  introduces  the  rainfall  model  equations,  and 
Subsection  2.2.3  gives  the  equations  for  the  channel  -  inflow 
and  rainfall  accumulator  states.  Section  2.3  presents  the 
discrete  equations  corresponding  to  Eq .  2.2-2. 

2.2.1  SSM  Model  Equations 

The  state  equations  corresponding  to  the  Sacramento 
Soil  Moisture  model  were  derived  from  the  NVSRFS  LAND  subrou¬ 
tine  documentation  (Ref.  8).  Six  states  (labeled  one  through 
six  in  Table  2.1-1)  are  necessary  to  represent  the  model:  two 
states  for  the  upper  zone,  three  for  the  lower  zone  and  one 
for  the  additional  impervious  area  content.  The  parameters  ot 
the  model  are  listed  in  Table  2.2-1.  The  relationship  given 
in  Table  2.2-1  between  the  instantaneous  drainage  coefficients, 
d^,  d^  and  d^'  and  their  counterparts  in  the  LAND  subroutine, 

UZK,  LZPK  and  LZSK  was  derived  in  Ref.  2. 

The  SSM  model  contains  several  t hreshol d- t ype  non- 
linearities  associated  with  the  bounds  on  the  contents  of  the 
elements  of  the  upper  and  lower  zones  and  with  constraints  on 


2-11 


THE  ANALYTIC  SCIENCES  CORPORATION 


TABLE  2.2-1 
SSM  MODEL  PARAMETERS 


T-  T'  it' 


NOIATION 

1  .AN  1 » 

F.QU  IVAIKNT 

IH.S(  R  1 1’TION 

o 

xi 

IIZTWH 

Upper-zone  t  ension-w.K  ei  rapacity  (mm) 

*1 

1 1 7.  KWM 

Upper-zone  free-w.it  er  capacity  (min) 

l./.TWM 

Lower -zone  t  ens  i  nn-w.it  <*i  capacity  (nun) 

X'» 

1.7  K  I’M 

I.ower-zonr  primary  free-w.iter  capacity  (mm) 

X 

1.7.  f  SM 

I  .owe  i  -  zone  secondary  frre-watot  capacity  (mm) 

-  f  n(  1  -UZK>/2<» 

Upper-zone  instantaneous  drainage  coefficient  (1  hi  ) 

-fn(  l-I.ZI'Kt/2'i 

Lower-zone  primary  instantaneous  drainage  coefficient  (I  /hi  ) 

■»; 

-(nl  1  -  t.ZSK  )/ZA 

l.owe  i -zone  secondary  instantaneous  drainage  coefficient  (1  hr) 

> 

7.1’KHl 

Parameter  in  percolation  function 

t* 

KKXI* 

Kxponettl  in  percolation  function 

"t 

PFRKF. 

Kt act  ion  of  percolated  waler  assigned  to  the  lower  zone 
free  water  aquifers 

p 

sim; 

fraction  of  base  Mow  not  appealing  in  liver  flow 

"i 

All  IMP 

fraction  of  basin  that  becomes  impervious  when  tension  water 

I equ i cement s  are  met 

■'2 

PCI  I N 

fraction  "1  basin  permanently  impervious 

f 

RSERV 

fraction  of  the  lower  zone  f fee-water  <  ,i|»,o  r  t  v  litravai  )  ab  i  e 
to  supply  lowei  zone  tension  i e«pi  i r ement s 

s 

R  I VA 

fraction  of  residual  evapnt  r  ansp  i  i  a  t  i  on  dem  mkI  a*  ln.il  lv 
contributed  hv  si  kmiii  sm  faces  and  iipaiiau  vegetation 

ratios  of  free  to  tension-water  content  for  the  two  zones .  A 
derivation  of  the  implications  of  these  const  lints  in  a 
continuous-t ime  SSM  model  was  included  in  Kef.  2.  For  the 
upper-zone,  the  constraint  can  be  written  as 


,  o  ^  ,  o 

x2/x2  -  xl/x] 


and  for  the  lower-zone. 


x4  +  x5 

1  -  — - 5  >  (1 

o  .  o  - 

XA  +  x5 


i  ■  1) 


(2.2-9) 


(2.2-10) 


Furthermore,  at  any  time,  t,  when  equality  holds  in  Kq .  2.2-9 
or  in  Eq .  2.2-10,  the  right  derivatives  of  the  states  must  satisfy 
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x2(t’ )/*2  -  xi ( c ' )/xi 


or 


x^(t+)  x,(t+)  +  X£.(t+) 

(1  -  r)  J  -  —  >  —  5 


x- 


o  .  o 

XA  x5 


(2.2-11 ) 


(2.2-12) 


respectively . 


In  addition  to  the  threshold- type  nonlinearities  there 
are  other  nonlinearities  in  the  SSM  model.  The  most  important 
nonlinearity  is  associated  with  the  percolation  function.  At 
any  given  time,  the  percolation  rate  from  the  upper  to  the 
lower  zone  is  given  by 


P 


P 


o 


1  + 


x3  +  x4  +  xj>\ 
o  ,  .o  .  o  ) 
x3  4  5/ 


with 


o  _ 
P  - 


X, 


+ 


d£ 


(2.2-13) 


( 2 . 2- 1 A ) 


Let  z^,  z^  and  z^  be  the  percolation  rate  inflow  into 
the  lower  zone's  tension-water,  primary  free-water  and  supple¬ 
mentary  free-water  elements  respectively.  In  the  SSM  model, 
under  normal  operation  (i.e.,  no  threshold  active),  t  he  total 
percolation  is  divided  into  the  lower  zone’s  elements  as 


Z3 


=  (1  -  Pf)P 


2/  = 


O  ,  O 

x4  x5 


20  ~ 

2  .  —  -  — i 

o  o 

XA  X5 


PfP 


(2.2-15) 

(2.2-16) 
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=  1  - 


o  ,  o 
XA  +  x5 


2d  -  x4/xj) 

2  -  —  -  ^ 
o  o 


PfP 


(2.2-17) 


Thus,  a  fraction,  (1  -  p^),  of  the  total  percolation  is  as¬ 
signed  to  the  tension-water  element  and  the  remainder  of  the 
percolation  is  divided  between  the  free-water  aquifers  accord¬ 
ing  to  Eqs.  2.2-16  and  2.2-17. 

In  Refs.  2  and  3  the  exact  differential  equations 
of  the  SSM  model  were  derived.  The  equations  given  there  cor¬ 
respond  to  the  distribution  of  the  percolation  rate  indicated 
by  Eqs.  2.2-15,  2.2-16  and  2.2-17.  There  is  a  difficulty  asso¬ 
ciated  with  the  use  of  this  distribution  of  the  percolation 
rate:  if  ,  the  fraction  of  p^p  assigned  to  the  primary 

free-water  aquifer,  Eq .  2.2-16,  is  larger  than  1  and  the  frac¬ 
tion  assigned  to  the  secondary  free-water  aquifer,  Eq .  2.2-17, 
is  negative  for  a  certain  range  of  values  of  x^  and  x^.  Ac¬ 
cordingly  the  distribution  of  the  percolation  rate  to  the  lower 
zone  was  slightly  modified.  The  modification  (Ref.  A)  is 
described  below. 


The  percolation  function,  Eq .  2.2-13  can  also  be 
written  as 


p  =  <4 


x?  —  +  d"  x?  —  +  p 
A  o  £  5  o  K 

X2  X2 


(2.2-18) 


with 


X3  *  XA  +  x5 \ 
o  ,  o  ,  o  I 
x3  XA  x5  / 


(  2  .  2  -  1  4  ) 


2- 1 A 


With  the  aid  of  Eq .  2.2-18  the  total  percolation  rate  can  be 
interpreted  as  consisting  of  two  parts.  One  part,  correspond¬ 
ing  to  the  first  two  terras  in  Eq .  2.2-18,  depends  on  the  maxi¬ 
mum  baseflow  rate  from  the  free  water  aquifers  and  on  the 
availability  of  water  in  the  upper-zone  free-water  element. 

The  other  part,  corresponding  to  the  term  p  in  Eq .  2.2-18, 
depends  on  the  lower  zone's  deficiency  ratio  as  well  as  on  the 
upper-zone  free-water  normalized  content. 


The  modified  distribution  of  the  total  percolation  is 

:  (1  -  Pf)  p  (2.2-20) 


=  d)  x. 


x2 

+ 

x4 

x4 

o 

x2 

o 

x4 

‘  x4 

+ 

o 

x5 

x5 

x2 

-f 

o 

x5 

- 

x5 

x2 

o 

x4 

-  x4 

+ 

o 

x5 

*  x5 

(2.2-21  ) 


(2.2-22) 


instead  of  Eqs.  2.2-15,  2.2-16  and  2.2-17.  Thus,  each  of  the 
free  water  aquifers  receives  a  part  of  the  percolation  which 
is  proportional  to  its  maximum  outflow  rate  and  to  the  availa¬ 
bility  of  water  in  the  upper-zone  free-water  element,  and  a 
second  part  which  depends  on  the  element's  deficiency  as  a 
fraction  of  the  total  deficiency  in  the  free  water  elements 
and,  also,  on  the  lower  zone's  deficiency  ratio. 

With  the  aid  of  the  functions 


hf(q)  = 


j  1  if  r)  >_  0 

I  o  i  f  q  <  o 


(2.2-23) 


he<n )  =  i  -  hf (n 


(2.2- 24  ) 
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The  state  equations  of  the  SSM  model  with  the  modified  distribu¬ 
tion  of  the  percolation  rate  to  the  lower  zone  can  be  written 
as  : 


*1  =  (“1  -  “2  |) 


—  )  lhe(-xi  '  x°>  +  hf(xi  "  Xi)  hf(u2  *  u]  )  1 


kp  - 

\x 2  xi  /  \x; 


hf(g* ) 


o  .  o 

X1  +  x2 


U1  "  u2  o 
X1 


—  -  d,  x,  -  p  h 


fl  o 

\x2 


|he  ( g  ’  ) 


(2.2-25) 


=  l(uj  -  u2)  hf(x1  -  Xj)  he(u2  -  Uj )  -  dux2  -  pi 
x  ^he^x2  ‘  x2^  +  hf^x2  "  x2  ^ 

*  hf  ( dux2  P  ~  ( u  ^  ”  u2)  h  j"(Xj-Xj  )  h^  (  u2  -  u  j  |  1 


Jh/^.!i\+hp.^\h 

\2  xl/  \x2  X1  / 


0 

i2 

0 

.  O 

X1 

+  x2 

ui  -  u9  —  -  d  x9  -  p  h 


(4  ‘  4)vs 


(2.2-26) 


-u2(l  -  xj/x°)  -T~~  +  U  -  Pf)P 

X1  x3 


X  |ht_(x2  -  x° )  +  hj-(x-i  -  x”)  hj-(m)  |  |h^(w)  *  h((v)h^(L" 


+  ‘u2(1  *  xl/xl)  +  * 

X1  X  3 
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(he(x3  -  x ^ )  +  hj(Xg  -  x° )  h f ( m ) ]  hf(w)  hf(g") 

(2.2-27) 


r  o 

X.  X.  -  X. 

x,  =  -d '  x ,  +  d'x°  —  +  - - - 

4  £4  £  4  o  o  ,o 

L  X2  X4  '  X4  X5  X5 


X  lhe(*3  *  x°)  +  hf(x3  -  x°)  hf(m))[he(w)  +  hf(w)  h^Cg")) 


r  o 

„  X  x.  -  X. 

+  -d’x.  +  d’x.  —  + - - - - — 

£4  £  4  o  o  ,o 

L  X2  X4  ‘  X4  +  X5  '  X5 


(p  -  z) 


x  Ihe(x3  *  x°)  +  hfCx3  "  *3)  hf(m)]  hf(w)  hf(g”)) 


/  0 

X  X.  -  X. 

>  H'r  t  H'v°  2  1  4  4 

p  -  u2(l  -  Xj/x,)  o  o 

)  £  4  £4o  0  0 

(  x2  X4  '  X4  +  X5  '  X5 

L  X1  x3. 

x  h.(x_)  h  (m) 
I  J  e 


(2.2-28) 


r  o 

X  x  -  X 

x  —  -d"x  +  d"x°  —  + - p  p 

5  £  5  £  5  o  o  o  f 

L  x2  X4  ~  X4  +  X5  X5 


X  lfie(x3  -  x“)  +  hf(x3  -  x3>  hf(m) )  (h^(w)  +  hf(w)  he  (g" ) ) 


r  o 

x0  x_  -  xc 

+  -d''x,.  +  d”x^  —  +  — — — - - - - 

£  5  £  5  o  0  .0 

L  X2  X4  '  X4  X5  ‘  X5 


( p  *  z  ) 


x  ^he^x3  "  x3^  +  hf*-x3  “  x3 lif(n»)]  hf(w)  hf(g") 


1  x„  x_  -  xc  x  .  ) 

/  £  5  £5o  o  o  2  1  1  t>  o  j 

<  x2  x4  -  x4  +  x5  ~  x5L  x  1  *  x 3  J  I 


X  h  (xj  h  (m) 
f  3  e 


(2.2-24) 
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- 


(u,  -  u0)  hr(x,  -  x?)  -  u0(l  -  X,/X?) 


j 2 1  1  f v  1  \ 

2 

-  <U1  •  U2){^J  hfUl  '  X°) 


l'-'l'  o  .  o 
X1  x3 


-  (uj  -  u2  -  dux2  -  p) 


XAi 

1  -1-5 

x3> 


h ^ ( x -  x  ^  )  h £ ( x 9 


o  x 
x9) 


(2.2-30) 

In  the  above  equations  and  u2  stand  for  the  instantaneous 
precipitation  and  potential  evapot ranspi rat  ion  demand  rates, 
respectively.  In  addition 

8>  =  ^o  ful  "  u2  ~oVhe(xl  "  xl)  +  hf(xl  '  xl)  hf(u2  U1  }  1 

X1  \  xl/ 

*  ^  t  <u1  -  u2)  hf(x1  -  x°)  h0(u2  -  Uj)  -  dux2  -  pi 


x  {h^(x9  -  x° )  +  hf(x2  -  x°) 


x  hf[dux9  +  p  -  (uj  -  u 2 )  hf(Xj  -  x^)he(u2  -  u ^ ) ] } 

(2.2-31  ) 


w  =  (1  -  r ) ( 1  -  x3/x°)  - 


1  - 


x4  *  x5 
o  .  o 
x4  +  x5 


(2.2-32) 


1  l  X5 

z  =  - — - — — -  ill  -  r)  (x.  +  x,)  u  (1  -  x  /x  )  - - 

o  .  , ,  w  o  ,  o.l  U  5  2  11  o  o 

x^  +  (1  -  r ) (x^  +  x^)’  x  +  x 


+  x3(p  ‘  d£x A  "  d£x5)( 


(2.2-33) 
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g"  =  z  -  (1  -  pf)p 

and 


(2.2-34) 


m  =  u2(l  -  Xj/x®)  Q  3 — -  -  (1  -  pf)p  (2.2-35) 

X1  +  x3 

The  propagation  of  state  values  in  the  simulation 
model  and  the  EKF  technique  are  both  based  in  successive  lin¬ 
ear  approximat ions  to  the  state-equations.  For  a  short  time 
interval  the  nonlinear  vector  function  F^  in  Eq .  2.2-1  is  re¬ 
placed  by  its  tangent  hyperplane  at  the  point  corresponding  to 
the  state  value  at  the  beginning  of  the  time  interval.  The 
choice  of  the  length  of  the  time  interval  is  discussed  in  Sec¬ 
tion  2.4.  The  inference  is  that  the  present  approach  to  simu¬ 
lation  and  filtering  imposes  restrictions  on  the  form  of  the 
state-equations:  the  nonlinear  function  £c  must  possess  con¬ 

tinuous  partial  derivatives  with  respect  to  all  states. 

The  function  F^  that  results  from  using  the  differen¬ 
tial  equations  2.2-25  through  2.2-30  directly  in  Eq .  2.2-1 
does  not  have  continuous  derivatives.  In  fact,  it  is  not  even 
continuous.  The  discontinuities  arise  from  the  threshold  val¬ 
ues  associated  with  the  elements  of  the  upper  and  lower  zones 
as  well  as  from  the  constraints  on  ratios  of  free  to  tension- 
water  present  in  the  SSM  model. 

Use  of  the  EKF  technique  on  the  reduced  SSM  model 
requires  modification  of  the  stat e-equat ions  2.2-25  through 
2.2-30  to  eliminate  the  threshold  discontinuities.  This  is  a 
step  of  cardinal  importance  in  the  analysis.  Some  of  the  most 
essential  features  of  the  SSM  model  such  as  the  supply  of  water 
to  the  upper-zone  fret-water  element,  the  distribution  ot  water 
to  the  lower  zone  and  the  appearance  of  surface  runoff  <it< 
critically  affected  by  the  approximation.  Arbitrary  smoothing 
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of  the  discontinuities  in  the  state-equations  may  yield  to 
physical  inconsistency  (the  principle  of  mass  conservation  may 
be  violated)  and  mathematical  inconsistency  (a  solution  may 
fail  to  exist  )  . 

An  analogy  with  electric  circuits  was  used  in  Ref.  5 
to  obtain  a  set  of  smoothed  equations  corresponding  to  the  SSM 
model  with  the  constraints  on  ratios  of  free  to  tension  water 
removed  from  the  model.  It  is  shown  in  the  Appendix  that  for 
all  basins  examined  the  constraint  for  the  upper  zone  is  super¬ 
fluous.  Even  though  situations  in  which  the  lower-zone  con¬ 
straint  is  activated  cannot  be  ruled  out  a  priori,  they  have 
not  been  observed  in  the  results  obtained  to  date.  When  the 
lower-zone  constraint  is  introduced,  the  smoothed  equations  of 
the  SSM  model  become 


moisture  evapot ranspi rat  ion  excess  upper-zone 
input  from  the  upper  zone  tension-water  supply 


(2.2-36) 


percolation  to  lower  zone 

(2.2-37  ) 
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X3  = 


( 1-Pf )P 


percol at  ion 
supply 


~5~-°  O-Xj/xJtxj 

■VT'T 


g(x3  ) 


evapotranspi rat i on  excess  lower-zone 
from  the  lower  zone  tension-water  supply 


+  f (z , ( 1  -  Pf)p)  h(w) 

effect  of  lower-zone 
constraint 


(2.2-38) 


o 

x4  ~x4 


I  PfP  +  g(>:3  ,x°)  ]  +  x 


o  ,  o 
x4’x4+x5”x5 


X  o 
o  2 

-'4  o 


moisture  supply  (percolation  +  excess  tension-water) 


X  r 


d;  x4 


baseflow  from 
primary 


o 

X5~X5 

x4-x4+x?-x5 


x4'x4 


o  ,  o 
x4"x4+x5~x5 


f (z  ,  ( 1  -  pf)p) 


h  ( w ) 


effect  of  lower-zone  constraint 

(2.2-39) 


x? 

PfP  +  g(x3,x°)]  +  d;’  X?  — 


‘S  -'5  .o 
‘0 


moisture  supply  (percolation  +  excess  tension-water 


di‘  x5 


baseflow  from 
supplementary 


o 

x  5  "  X  5 


o  .  o 

x4*x4  x5*x5 


f  ( z  ,  (  1  -  p  f  )  p  )  h  ( w  ) 


effect  of  lower-zone  constraint 

(2.2-40) 


where  w  and  z  are  given 
and  where  the  functions 


by  Eqs .  2.2-32 
f.  g  and  h  are 


and  2.2-33  respectively 
given  by 
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0 

if 

n 

<  rj  -6 

f  ( o .  n  >  =  \ 

(n+6-fj  )2 
46 

if 

n  ■ 

-6  <  rj  <  rj+6 

I 

n-n 

if 

n 

>  rj  +6 

g(n.n°)  =1 

i  ° 

i  e  .  o 

(n*n 

>2 

if 

if 

^  o 

n  £  n 

v  o 

n  >  n 

ln 


h(n> 


0 

(n  +  <5/2  ?2 
6 

in  -  fi/2)2 

6 

1 


if  i)  <  -6 
if  -6  <  n 

if  -6/2  <  n 
i  f  n  1  0 


<  -6/2 

<  0 


(2 


(2 


(2 


where  6  and  e  are  constants,  typical  values  of  which  are 
=  0.01,  e  =  100.  Graphs  of  the  functions  f,  g  and  h  are 
in  Figs.  2.2-1,  2.2-2  and  2.2-3,  respectively. 


R-51647 


2-01 ) 


2-02) 


2-03) 


6 

given 


F igure  2.2-1 


Graph  of  the  Function  f(r|«n) 
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bounds 

ceeded 


R-51646 


Figure  2.2-2 


Graph  of  the  Function  g(n.H°) 


R -595  16 


h  ( rj  ) 


Figure  2.2-3  Graph  of  the  Function  h(q) 


As  a  consequence  of  the  approx i ma t i on  the  upper 
on  the  state  values  niaj7 ,  on  occasion,  be  slightly  ex 
Examples  of  this  behavior  are  given  in  Section  2.6 
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By  increasing  the  value  of  the  constant  e  in  Eq  .  2.2-42,  the 
possible  excess  of  state  values  over  their  bounds  may  be  re¬ 
duced.  Even  though  in  the  results  obtained  to  date  the  lower- 
zone  deficiency  ratio,  p,  defined  by 


p  -  1 


x3  ♦ 


+ 


o 

x5 


(2 .2-44 ) 


has  always  been  greater  than  zero,  there  is  the  possibility 
that  it  may  become  negative  during  extreme  flood  periods. 
This  ratio  appears  in  the  percolation  function  and  if  its 
value  is  negative,  p  given  by  Eq .  2.2-19  is  undefined.  The 
definition  of  p  is  completed  by  setting 


P 


p  >  0 


p  <  0 


(2.2-45) 


with  p  as  in  Eq .  2.2-44.  Thus,  when  the  lower  zone  deficiency 
ratio  is  negative,  water  is  drawn  out  of  the  lower  zone  into 
the  upper-zone  free-water  element.  If  the  latter  is  full,  the 
amount  of  water  drawn  out  of  the  lower  zone  contributes  to  the 
surface  runoff. 


The  equation  for  x^ ,  the  excess  of  the  additional  im¬ 
pervious  storage  over  the  upper-zone  tension-water  content,  is 


=  I g(Xj  ,Xj  )  -  g(x2 ,x£ ) )  1 1 


,  ,  o  .  2 

(xb/x3) 


u2  (  1 


X, 

/  O  \  V 

W  — ;--s 

xi  x:i 


(2.2-44) 
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The  equations  for  the  moisture  input  u^  in  Eq .  2.2-36 
are  discussed  in  the  following  section. 

2.2.2  Rainfall  Model 

The  Kalman  filtering  formulation  assumes  that  random 
inputs  driving  a  system  are  solutions  of  stochastic  differen¬ 
tial  equations.  The  simplest  such  input  is  the  so-called  first- 
order  Markov  process  (Ref.  7)  which  is  generated  by  the  equation 

Xg(t)  =  -  p  xg(t)  +  C(t)  (2.2-47) 

where  p  is  a  real  number  and  £  is  a  Gaussian  white  noise  proc¬ 
ess  with  constant  spectral  density  q.  If  Xg  is  interpreted  as 
a  rainfall  rate  with  units  of  mm/hr,  the  units  of  £  are  mm/hr^ , 
those  of  q  are  (mm/hr  )  /(1/hr),  and  the  units  of  p  are  1/hr. 

If  p>0  a  stationary  solution,  Xg(t),  to  Eq .  2.2-47  exists,  and 
its  covariance  function  is 

E [ Xg ( t+x )  Xg(t) ]  =  e*PlTi  (2.2-48) 

The  quant  ity  1/p  is  called  the  correlation _ t  ime  of  the  process. 

The  process  Xg(t)  has  a  mean  value  of  zero.  Its  sam¬ 
ple  functions  oscillate  about  the  time  axis,  and  even  though 

1  /2 

values  much  larger  than  the  standard  deviation,  (q/2p)  , 

have  a  low  probability  of  occurrence,  there  is  no  actual  upper 
or  lower  bound  on  the  values  the  process  may  take.  A  precipi¬ 
tation  model  should  reflect  the  existence  of  periods  in  which 
there  is  no  rain  and  the  fact  that  rainfall  rates  are  never 
negative.  Thus,  Eq .  2.2-47,  per  se ,  does  not  constitute  an 
appropriate  precipitation  model  in  a  simulation.  However,  a 


*  E  denotes  the  mathematical  expectation  operator. 
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simple  nonlinear  t rans format  ion  of  the  process  Xg(t)  yields  a 
useful  precipitation  model  for  simulation  purposes.  Consider 
the  precipitation  rate  defined  as 

Uj ( t )  =  f(x8(t),xg)  (2.2-49) 

where  the  function  f  is  given  by 

l  0  i  f  r)  £  H 

f (n  .n  )  =  (2.2-50) 

(  n-n  if  n  >  n 

Whenever  the  values  of  the  process  Xg  are  less  than  the  thresh 
old  value  Xg,  the  precipitation  rate  is  zero.  If  Xy  is  larger 
than  Xg,  their  difference  is  taken  as  the  precipitation  rate. 
An  example  of  a  typical  realization  of  six-hour  accumulated 
rainfall  obtained  with  this  model  is  given  in  Fig.  2.2-4. 

Other  examples  are  given  in  Section  2.6. 


Six-Hour  Accumulated  Precipitation 


Figure  2.2-4 


Example  of  Accumulated  Precipitation  Produced 
with  the  Rainfall  Simulation  Model 
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The  fraction  of  time  during  which  rain  occurs  is  equal 
to  the  probability  that  Xg(t)  is  greater  than  Xg  and  is  given  by 


Pluj(t)  >  0]  =  1  -  <)»(Xg/o) 


(2.2-51  ) 


where 


♦  (a)  =  4=  f 

y/2n  -OB 


a  2  ,,, 

-s  /2  . 
e  ds 


(2.2-52) 


is  the  standardized  normal  distribution  function,  and 


a  -  (q /2p) 


(2.2-53) 


is  the  standard  deviation  of  Xg(t).  The  expectation  and  vari¬ 
ance  of  the  precipitation  rate  given  that  it  is  actually  rain¬ 
ing  can  be  computed  to  be 


E| Uj ( t )  Uj ( t  )  >  0 |  =  Ao  -  Xg 


u^ ( t )  Uj ( t )  >  0 


where 


=  a2 [ 1-A2 |  +  Aoxg 


(2.2-54) 

(2.2-55) 


-2  /0  2 
-Xq/2o 

A  =  e  ° _ /n  n 

1  -  o ( x  g/n  ) 


(2 .2-56) 


Equations  2.2-51,  2.2-54,  and  2.2-55  may  be  used  in 
fitting  the  model  parameters,  p,  q,  and  Xg ,  to  data  records. 
The  seasonal  variation  of  rainfall  can  be  modeled  by  letting 
the  model  parameters  vary  with  time.  For  the  results  included 
in  Section  2.6  corresponding  to  the  Bird  Creek  basin  and  for 

the  precipitation  record  of  Fig.  2.2-4  the  correlation  time, 

v  v 

1/p,  was  taken  to  be  one  day,  q  -  0.285  ( min/hr4"  )  /  (  1 /'hr )  and 
xk  -  2  mm/hr.  And  for  the  results  corresponding  to  the  White 
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River  basin  in  Section  2.6,  the  parameters  are  1/p  =  12  hr,  q 
=  0.015  (mm/hr^ ) ^/( 1/hr )  and  Xg  =  0.153  mm/hr.  These  parame¬ 
ters  were  chosen  on  the  basis  of  NWS  values  of  mean  areal  pre¬ 
cipitation  (MAP)  for  the  Bird  Creek  and  White  River  basins 
corresponding  to  the  periods  between  April  19,  1959  and  May 
19,  1959  and  between  May  1,  1968  and  June  1,  1968,  respectively. 
These  periods  were  chosen  because  the  relatively  high  quanti¬ 
ties  of  prec ipi tat  ion  that  were  recorded  served  to  illustrate 
the  nonlinear  behavior  of  the  model  in  the  results  presented 
in  Section  2.6. 

The  function  f  in  Eq .  2. 2-49  does  not  have  a  deriva¬ 
tive  when  Xg  =  Xg.  Since  the  operation  of  the  simulation  model 
requires  the  existence  of  derivatives  for  the  linearization  of 
the  equations,  the  function  f  was  replaced  by  f  as  defined  by 
Eq.  2.2-41. 


The  above  model  is  useful  only  for  simulation  pur¬ 
poses.  It  is  referred  as  rainfall  simulation  model  in  the 
sequel.  For  reasons  that  were  discussed  in  Ref.  5,  this  model 
is  inappropriate  for  the  filtering  formulation. 

Instead  of  the  above  simulation  model,  the  model  used 
in  the  filtering  formulation  is  a  linear  rainfall  rate  model 
that  has  the  same  first  and  second-order  moments  as  the  rainfall 
simulation  model.  The  rainfall  rate  is 

u j ( t )  =  u j  +  Xg( t  )  (2.2-57) 

where  U|  is  the  average  rainfall  rate  given  by 

Uj  =  P | U] ( t )  >  0|  E|Uj(t)ju](t)  >  01  (2.2-58) 
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or,  in  terms  of  the  quantities  previously  introduced,  by 

uj  =  [1  -  0<Xg/a>]  |Ao  -  x8J  (2.2-59) 

The  process  x^(t)  is  a  first-order  Markov  process  satisfying 
the  stochastic  differential  equation 

Xg  ( t )  =  -  p  Xg  ( t )  +  C'U  )  (2.2-60) 

where  £'(t)  is  a  Gaussian  white  noise  process  with  spectral 
density 

q'  =  j|l  -  0  ( Xg/ CT  )  ]  (1  +  (Xg/o)2  -  A.\g/o  ] 

-  (1  -  0(Xg/o)]2  [A  -  Xg/ol2|q  (2.2-61) 

The  values  of  u^  and  q'  associated  with  the  rainfall 
simulation  models  previously  introduced  are  u^  =  0.132  mm/hr 
q'  =  0.0165  ( mm/hr2 ) 2/( 1/hr )  for  the  Bird  Creek  basin  and  Uj  = 
0.0583  mm/hr,  q'  =  0.00259  (mm/hr2 )2/(  1  /hr )  for  the  White  Rivet- 
basin.  The  rainfall  model  corresponding  t o  Eqs.  2.2-57  and 
2.2-60  is  called  rainfall  filtering  model  in  the  sequel. 

2.2.3  Equations  for  the  Channel  -  1 n  f 1 o  w  and 
Rainfall  Accumulator  States 


The  contributors  to  the  channel  inflow  accumulator 
are:  rainfall  over  permanently  impervious  areas,  direct  run¬ 

off  from  the  additional  impervious  area,  surface  runoff  from 
the  upper  zone  and  from  the  additional  impervious  area,  inter¬ 
flow  and  base  flow.  To  all  these  contributions  t he  evapot rans- 
pi rat  ion  from  stream  surfaces  and  riparian  vegetation  is  sub¬ 
tracted.  The  equation  for  x^  is  easily  found  to  be 

x?  =  f(y,0)  (2.2-62) 
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where  f  is  given  by  Eq .  2.2-41  and 


21  1,KV  1’  1 


+  (l-aj-a2) 


d  x0  + 


+  (l-a2>g(x2 ,x^ ) 


-  s  u  2  <  1  -  x  i  /x 


X?)  (l  -  \ 

1  \ 


(2.2-43  ) 


The  equation  for  the  rainfall  accumulator  is,  simply. 

Xg  =  Uj  (2. 2  -(>4) 


2.3  MODEL  EQUATIONS  -  DISCRETE  TIME  COMPONENTS 

In  contrast  with  the  equations  lor  the  continuous  part 
of  the  model,  those  of  the  discrete  part  are  linear.  Let  the 
state-space  model  of  the  unit  hydrograph  be  (see  Section  3.4.4) 

*h(j)  =  %  -h(j"1)  +  Gh  uh(j)  (2.3-1) 

vh(j)  =  Hh  xh(j)  (2.3-2) 

where  x^  is  the  state  vector,  the  input,  v^  the  discharge 
rate,  <J>^  the  M*M  one-step  transition  matrix,  the  M*  1  input 
matrix  and  the  1><M  output  matrix  of  the  unit  hydrograph  and 
where  time  (j=0,l,...)  is  measured  in  multiples  of  tne  unit 
hydrograph  rate. 

It  is  convenient  to  express  the  discharge  at  time  j 
in  terms  of  the  state  vector  at  time  j-1.  Replacing  x^(j)  in 
Eq .  2.3-2  by  the  expression  on  the  right  of  Eq .  2.3-1  yields 
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vh(j)  =  Hh  +  Hh  Gh  uh(j)  (2.3-3) 

The  matrices  A  and  B  of  Eq .  2.2-2  can  be  easily  de¬ 
rived  from  Eqs.  2.3-1  and  2.3-3  as  follows.  Let  vk  be  a  crit¬ 
ical  time.  The  first  M  states  of  x^  ,  those  corresponding  to 
the  unit  hydrograph,  are  governed  by  Eq.  2.3-1  when  j  is  re¬ 
placed  by  (vk)°,  j-1  by  (vk)  ,  and  u^(j)  by  Xy|(vk)  ),  the 
latter  being  the  current  value  of  the  accumulated  channel  in¬ 
flow.  The  last  state  of  x^  ,  the  accumulated  channel  discharge 
at  time  (vk)  °,  can  be  written  as  the  sum  of  the  previous  con¬ 
tents  of  the  accumulator,  x^Q+^[(vk)  J,  and  the  present  con¬ 
tribution  from  the  channel  routing  system  given  by  Eq .  2.3-3, 
i  •  e  .  , 


X10+M[(vk)°]  =  X10+M[(vk)’1  +  Hh  *h  ^h1(vk) 


+  Hh  Gh  x7((vk)'J 


(2.3-4) 


The  matrices  A  and  B  are 

— M - 1 

f 

I 

A  = 


\-l-\ - 2—1 


B  = 


0  . 

•  •  0 

:  o 

0\ 

■ 

Gh  ' : 

• 

0  • 

•  •  6 

:  o 

0 

\o  . 

•  •  0 

IU.  o 
h  h> 

°/ 

I 

M 


(2.3-5) 


(2.3-6) 


2-31 


THE  ANALYTIC  SCIENCES  CORPORATION 


2.4  OPERATION  OF  THE  MODEL  IN  THE  SIMULATION  MODE 

Between  two  consecutive  critical  times  (v-l)k  and  vk. 
The  equations  for  the  continuous  part  of  the  system, 

xc(t)  =  Fc(xc(t),t>  +  GcC(t)  (2.4-1) 

are  linearized  at  the  operating  point  corresponding  to  time 
tQ=vk.  The  resulting  linear  equations 

x  ( t )  =  F  ( x  ( t )  -  x  ( t  ) 1  +  h  +  G  £ ( t )  (2.4-2) 

— c  c  —  C  —C  O  -  C 


wi  th 


F 

c 


3  F 


3x 

-c 


~  ( X  (  t  )  ,  t  ) 

T  -c  o  o 


(2.4-3) 


and 


h  =  F  (x  (t  ) ,t  )  (2.4-4) 

-  -c-c  o  o 

are  used  to  obtain  a  first  approximation  to  the  state  values 

at  time  t  +A  with  A  =  vk  -  (v-l)k  =  k. 
o 

The  state  vector  at  the  end  of  the  time  interval  of 
length  A  is  the  value  of  the  solution  to  Eq .  2.4-2  at  time 
1 0+A .  This  value,  denoted  by  x^(to+A),  is  found  to  be 

-c(t0+A)  =  ^c(to)  +  b  +  £'  (2.4-3) 

where  I  is  the  identity  matrix,  <t>(,  (s)  the  transition  matrix 
for  an  interval  of  length  s  given  by 


Vs> 


F  s 
c 


(  2 . 4  -  (. ) 
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and  ^ '  a  vector  randomly  chosen  from  the  multivariate  normal 
distribution  with  mean  zero  and  covariance  matrix 

A 

Q  =  q  /  ♦  (A-s)  C  G^(i-s)  ds  (2.4-7) 

~0  c  c  c  c 

The  vector  '  accounts  for  the  contribution  of  the  random  input  , 
£  ,  of  the  rainfall  simulation  model  to  the  change  in  the  state 
values  during  the  time  interval  of  length  A. 

Because  of  the  strong  nonlinearities  present  in  the 
model,  it  is  possible  that  the  behavior  of  the  linearized  sys¬ 
tem,  Eq .  2.4-2,  departs  considerably  from  that  of  » he  nonlinear 
system,  Eq .  2.4-1,  during  the  time  of  propagation;  i.e.,  the 
length  of  the  time  interval.  A,  may  be  inappropriately  large 
for  a  valid  approximation.  In  order  to  avoid  this  difficulty, 
the  length  of  the  time  interval  for  propagation  of  the  linear 
system  is  chosen  adaptively.  The  simulation  utilizes  a  nominal 
step  length  of  A  =  k  hr,  but  if  any  of  a  collection  of  inequal¬ 
ities,  W(x.c».x^)  £  0,  is  violated,  the  length  of  the  time  inter¬ 
val  is  halved  and  propagation  is  attempted  anew.  When  propa¬ 
gation  for  a  subinterval  is  successfully  completed,  Eq .  2.4-1 
is  linearized  at  the  new  operating  point,  A  is  set  equal  to 
the  time  step  necessary  to  reach  the  next  critical  time  and 
the  procedure  is  repeated. 

Four  inequalities  were  used  in  the  simulation.  The 
first  two  place  limits  on  the  changes  in  the  contents  of  the 
upper-zone  tension-water  and  free-water  elements;  these  dit- 
ferences  are  not  allowed  to  exceed  1  mm  ;  i.e. , 

|*i“„>  •  *i<vi>  I  "  1 

*This  amount  is  a  conservative  choice.  Good  results  were 
obtained  with  its  use. 
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|X2 ( tQ )  ~  X^(tQ+A)  |  <  1  (2.4-9) 

The  third  inequality  concerns  the  percolation.  It  is  required 
that 


P ( x ( t  ) )  +  p(x ' ( t  +A  )  ) 

- 2 ° -  a  <  1  (2.4-10) 


The  expression  on  the  left  in  Inequality  2.4-10  is  Simpson's 
integration  formula  applied  to  the  instantaneous  percolation 
rate,  p.  Thus,  2.4-10  can  be  interpreted  as  a  limit  on  the 
total  percolation  during  the  time  interval  of  length  A.  The 
last  inequality  sets  a  limit  on  the  amount  of  surface  runoff: 


[g(x2( tQ)  ,x°)  +  g(x^(tQ+A) ,x°)  ] 


A  <  1 


(2.4-11 ) 


When  the  critical  time  vk  is  reached,  the  last  value 
of  the  continuous  state  vector,  x  ,  obtained  by  the  propagation 
of  Eq .  2.4-1  is  interpreted  as  x^Kvk)  ).  The  corresponding 
value  of  the  discrete  state  is  x^|(vk)  |  =  x^ [ ( ( v -  1 ) k  )  r  |  .  The 
discrete  part  of  the  model  is  updated  next  according  to 

xd!(vk)°]  =  A  xd((vk)"j  +  B  xr|(vk>-]  (2.4-12) 

where  the  matrices  A  and  B  are  given  by  Eqs.  2.3-9  and  2.3-6. 

At  this  point  in  the  simulation  the  output  measurements, 
Eqs.  2.2-7  and  2.2-8,  are  evaluated.  It  is  convenient  to  intro¬ 
duce  a  matrix  notation  for  this  computation.  Equations  2.2-7 
and  2.2-8  can  be  represented  as 


y(vk)  =  H(vk)  x((vk)°|  +  £  ( v  k  ) 


(2.4-13) 
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where  the  quantities  y,  H  and  4  are  given  by 


when  v  is  a  multiple  of  m £;  i.e.,  when  there  are  rainfall  and 
discharge  measurements,  and  given  by 


4  =  (2.4-19) 

when  v  is  not  a  multiple  of  but  is  a  multiple  of  £  ,  i.e., 
when  there  is  only  a  rainfall  measurement.  If  v  is  not  a  mul¬ 
tiple  of  £  there  are  no  measurements  computed.  In  the  above 
equations  and  ^  represent  measurement  noise.  Their  values 

are  independent  samples  from  the  normal  distribution  with  zero 

2  2 

mean  and  variances  and  o 2-  respectively. 

Since  the  filter  is  not  in  operation,  the  definition 
of  time  (vk)  is  superf1  ious  in  the  simulation  mode.  Equiva¬ 
lently  , 

x( (vk)+  |  =  x | ( v  k  )  °  |  (2.4-20) 
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As  a  final  step  of  the  state  propagation  correspond¬ 
ing  to  a  critical  time,  the  contents  of  the  appropriate  ac- 


cumulators  are 

reset  to 

zero ,  i . e . , 

x7[ (vk)rl 

=  0 

;  v  =  0 ,  1  , 

2,  ... 

(2.4-21  ) 

x<j[  (vk)r] 

=  0 

O 

II 

> 

2£,.  .  . 

(2.4-22) 

X10+Ml(vk) 

II 

O 

< 

If 

o 

3 

Kb 

,  2m£ .... 

(2.4-23) 

and  propagation  of  the  continuous  part  of  the  state  vector  to 
the  next  critical  time  begins. 

2.5  OPERATION  OF  THE  MODEL  IN  THE  FILTERING  MODE 

The  operation  of  the  filter  follows  in  the  same  lines 
as  the  simulation  procedure  described  in  Section  2. A.  The 
system  equations  are  identical  to  the  equations  used  in  the 
simulation  mode  except  that  the  rainfall  filtering  model  is 
used  instead  of  the  rainfall  simulation  model.  The  state- 
estimates,  x,  are  propagated  according  to  the  equations 

f£c(t>  =  Fc(xc(t),t)  ;te((v-l)k,vk)  (2.5-1) 

xd(t)  =  0  ;te((v-l)k,vk)  (2.5-2) 

xdl(vk)°J  =  A  x d ( ( v  k  )  I  +Bxc((vk)~]  (2.5-3) 

However,  in  order  to  compute  the  Kalman  gains,  it  is  necessary 
to  propagate  the  state  error  covariance  matrix,  P(t),  defined 

by 


P(t)  =  E{ [x(t  )  -  x( t )) 1  | x ( t )  -  x( t ) )T} 


(2.5-4) 
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Following  the  subdivision  of  the  state-vector  into 
its  continuous  and  discrete  parts,  xc  and  xd  ,  define 


Pc(t)  =  E{(xc(t) 

-  xc(t)] 

(xc(t)  - 

xc(t)]T! 

fcd(t)  -  E{(xc(t) 
Pdc<l)  =  pId(t> 

-  xc ( t ) J 

(xd(t)  - 

xd(t)]Ti 

Pd(t)  =  E{ (xd(t) 

-  xd(U] 

(xd(t)  - 

xd(t)]Ti 

Thus  P(t)  is  partitioned  as 


P(t) 


( 


pdc(t)  :  pd(t>  / 


(2.5-5) 


(2.5-6) 

(2.5-7) 

(2.5-8) 


(2.5-9) 


The  operation  of  the  filter  is  best  described  in  the 
form  of  an  algorithm. 


Filtering  Algorithm 

1.  (Initialize  elapsed  time,  state-estimates  and  error  covari¬ 
ance  matrix  | 

Set  t  «-  0,  x  +  x(0),  Pc  -  Pc(0),  Pcd  -  Pcd(0),  Pd  -  P()(0) 

2.  (Initialize  step  size.  A,  and  time  to  reach  next  critical 
time,  A . | 

Set  A  <-  k,  A  <-  k 

3.  (Linearize  Eq .  2.5-1  about  current  operating  point  | 

3  F 

Compute  F  =  — £  (x ,  t) ,  h  =  F  (x  t) 

Bx 

-c 

A.  (Evaluate  transition  matrixl 

F  A 

Compute  4>  =  e 
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5.  [Compute  tentative  state-estimate] 

Set  x'  «-  x  +  (4>  -I)  F"1  h 
— c  — c  c  c  — 

6.  [Verify  suitability  of  step  size] 

If  none  of  the  inequalities  W(xc,x^)  <  0  is  violated  go  to 
step  8  (W  is  defined  by  Eqs.  2.4-8  through  2.4-11) 

7.  [Redefine  step  size) 

Set  A  *■  A/2  and  go  to  step  4 

8.  [Update  state-estimate,  elapsed  time  and  time  to  next 
critical  time ] 


9. 


Set  x  «■  x '  ,  t  «-  t+A,  A  ^  A  -  A 
-c  -c 

[Propagate  covariance  matrix  to  account  for  continuous 
transition ] 


Set  P  «•  <t>  P  4>  +  Q,  P  , 

c  c  c  c  cd 


<t>  P  , 
c  cd 


where  Q  is  the  matrix  defined  by  the  expression  on  the  right 
of  Eq .  2.4-7. 


10.  [Verify  if  critical  time  has  been  reached) 


If  A  /  0  set  A  «-  A  and  go  to  step  3 

11.  iPropagate  discrete  part  of  state  vector) 

Compute  state  vector  at  time  t°  by  setting 

x  ,  <-  A  x  ,  +  B  x 
-d  -d  -c 

12.  IPropagate  covariance  matrix  to  account  for  discrete 
t  rans i t ion  J 

Set  P  «-  A  P  .  AT  +  B  P  ,  AT  +  (BP.  AT)T  +  B  P  B1 
d  d  cd  cd  c 

and  P  ,  P  ,  A1  +  P  BT 
cd  cd  c 

in  the  above  order 


13.  [Determine  if  there  is  a  precipitation  measurement  | 
If  t  is  not  a  multiple  of  £k  then 
13.1  j Reset  channel  inflow  accumulator) 

State  at  time  tr  obtained  by  setting  Xy  •  0 
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13.2  [Modify  covariance  matrix  to  reflect  rebet) 

Set  the  seventh  row  and  column  of  P  equal  to  zero 

13.3  Go  to  step  2 

14.  (Determine  if  there  is  a  discharge  measurement] 

1ft  is  a  multiple  of  ml!  k  set  FLAG  *■  TRUE 

If  t  is  not  a  multiple  of  mi  k  set  FLAG  *■  FALSE 

(FLAG  =  TRUE  indicates  the  presence  of  a  discharge  measurement  ) 

15.  [Define  measurement  equations] 

If  FLAG  =  TRUE  define  H  as  in  Eq .  2.4-15  and  set 


If  FLAG  =  FALSE  define  H  as  in  Eq .  2.4-18  and  set 

x  +  yi  R  *  ai 

(yj  and  y2  are  the  values  of  the  rainfall  and  mean  discharge 

measurements,  R  is  the  covariance  matrix  of  the  errors  in 
the  measurements) 

16.  [Compute  Kalman  gains] 

Evaluate  K  =  P  HT(HPHT  +  R ) " 1 

17.  [Perform  Kalman  update  on  state-estimate! 

State  at  time  t  is  obtained  by  setting 

x  <-  x  +  K I  y  -  Hx  ] 

18.  [Modify  covariance  matrix  to  account  for  Kalman  update  | 

Set  P  «■  P  -  KHP 

19.  [Reset  accumulators] 

19.1  Set  x  —j  *  0  and  x(^  *  0 

19.2  If  FLAG  =  TRUE  set  xinjtM  •  0 

1  U  +  M 

2-39 


M 


THE  ANALYTIC  SCIENCES  CORPORATION 


20.  [Modify  covariance  matrix  to  account  for  reset) 

20.1  Set  the  seventh  and  ninth  rows  and  columns  of  P 
equal  to  zero 

20.2  If  FLAG  =  TRUE  set  the  last  row  and  column  of  P 
equal  to  zero 

21.  (Continue  operation] 

Go  to  step  2 

2.6  RESULTS 

The  parameters  of  the  SSM  model  for  the  Bird  Creek 
and  White  River  basins  were  used  in  testing  the  behavior  of 
the  simulation  and  filtering  algorithms  described  in  Sections 
2 . A  and  2.5. 

In  all  instances  data  were  generated  using  the  model 
in  the  simulation  mode.  These  data  were  then  used  by  the  fil¬ 
ter  to  obtain  six-hour  lead  forecasts  under  several  assumptions 
on  the  filter  parameters.  The  results  so  obtained  were  then 
compared  to  the  true  values  obtained  in  the  simulation.  Excel¬ 
lent  agreement  between  true  and  predicted  discharge  was  obtained 
in  all  test  cases . 

The  Bird  Creek  basin  model  was  used  to  investigate 
the  sensitivity  of  the  filter  behavior  to  changes  in  the  meas¬ 
urement  error  covariance  matrix,  R,  and  to  changes  in  the  ini¬ 
tial  state  error  covariance  matrix,  P(0).  The  White  River 
basin  model  was  used  to  illustrate  the  filter  response  for  a 
different  set  of  SSM  parameters. 

The  parameters  of  the  SSM  model  for  the  Bird  Creek 
and  White  River  basin  are  listed  in  Table  2.6-1.  In  addition 
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TABLE  2.6-1 

SSM  MODEL  PARAMETERS  FOR  BIRD  CREEK 
AND  WHITE  RIVER  BASINS 


PARAMETER 

BIRD  CREEK 

WHITE  RIVER 

.x  ^  ( mm ) 

120 

50 

x^  ( mm ) 

15 

30 

x°  ( mm ) 

160 

250 

x^  (mm) 

140 

80 

x°  ( mm ) 

14 

170 

du  (1/hr) 

1 .486E-2 

1 . 486E-2 

d^  (1/hr) 

5.452E-4 

5.875E-4 

d'1  (1/hr) 

5.612E-3 

2.578E-3 

y 

48 

25 

a 

2.1 

4 . 0 

Pf 

0.02 

0.4 

M 

3.55 

0.0 

al 

0.17 

0.0 

3  ^ 

0.001 

0.0 

r 

0.3 

0 . 3 

s 

0.0 

0.2 

to  the  parameters  given  in  the  table,  constant  rates  ol  u 
1.375E-2  mm/hr  and  =  4.583E-2  nun/hr  were  used  for  the 
stantaneous  potential  evapot ranspi rat  ion  demand  for  t he  B 
Creek  and  White  River  basins,  respectively. 
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Fifth-order  state-space  models  for  the  six-hour  unit 
hydrographs  for  the  Bird  Creek  and  White  River  basins  were 
used.  These  models  were  obtained  with  the  subroutine  REDO-LTHG 
(Ref.  1)  delivered  to  NOAA/NWS  for  use  as  a  constituent  of 
the  Version  5.0  NWSRFS  Forecast  Component.  The  theoretical 
basis  underlying  the  operation  of  this  subroutine  is  described 
in  Chapter  3  of  this  report.  Figure  2.6-1  compares  the  origi¬ 
nal  Bird  Creek  hydrograph  (solid  line)  with  the  approximation 
used  in  the  simulation  and  filtering  algorithms  (dashed  line). 
The  White  River  basin  hydrograph  (Fig.  2.6-2)  contains  five 
lags  and,  thus,  agrees  exactly  with  the  fifth-order  state- 
space  model  used  in  the  simulation  and  filtering  algorithms. 

For  the  Bird  Creek  basin  the  initial  state  values 
were  chosen  arbitrarily  as  Xj(0)  =  100  mm,  X2<0)  =  12  mm ,  x^(0) 
=  130  mm,  x^(0)  =  110  mm,  x^(0)  =  11  mm,  and  x^(0)  =  0  for 
i  >  5.  These  values  were  the  same  in  the  simulation  and  fil¬ 
tering  operations. 

For  the  Bird  Creek  model,  different  random  number 
sequences  were  used  in  testing  the  sensitivity  of  the  filter 
response  to  changes  in  the  measurement  and  initial  state  error 
covariance  matrices.  The  particular  precipitation  records 
included  in  this  section  are  not  typical  of  the  results  pro¬ 
duced  by  the  rainfall  simulation  model.  They  were  chosen  be¬ 
cause  the  substantial  amounts  of  rainfall  obtained  with  the 
particular  random  number  sequences  excite  extreme  dynamical 
responses  in  the  system  and  serve  to  illustrate  the  behavior 
of  the  filter  in  the  nonlinear  operating  region  for  the  upper- 
zone  elements. 

The  precipitation  record  of  Fig.  2.6-3  for  the  Bird 
Creek  basin  was  used  in  determining  differences  in  filter  re¬ 
sponse  to  changes  in  the  measurement  error  covariance  matrix, 

R.  This  prec i pi t a t i on  record  is  referred  as  BC1  in  the  sequel. 
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Bird  Creek  Unit  Hydrographs 


0  2  4  C  P 


TIME  (days) 

Figure  2.6-1  Six-Hour  Unit  Hydrographs  for  Bird  Creek: 

Original  Hydrograph  (Solid  Line)  and  Fifth- 
Order  State-Space  Approximation  (Dashed  Lint 


White  River  Unit  Hydrograph 


o  1 


T I  Ml  (days) 

Six-Hour  Unit  Hydrograph  for  White  River  Basin 


F igur e  2 . 6-2 
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Six-Hour  Accumulated  Precipitation 
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Figure  2.6-3  Precipitation  Record  BC1 

For  the  comparisons  of  filter  performance  under  dif¬ 
ferent  assumptions  on  the  measurement  noise  level,  simulated 
values  of  accumulated  prec i pi t a t i on  and  i ns t an t aneous  discharge 
were  produced  at  six-hour  intervals.  Thus,  in  the  notation  of 
previous  sections,  the  values  of  k,  £,  and  m  were  k=6,  £=1  and 
m=l  . 


Results  for  two  different  measurement  error  covari¬ 
ance  matrices  are  illustrated  in  Figs.  2.  b- A  through  2.6-10. 
Solid  lines  represent  the  values  obtained  in  the  simulation 
while  dashed  lines  correspond  to  the  one-step  ahead  (six  hour) 
filter  prediction.  Even  numbered  figures  correspond  to  the 
measurement  covariance  matrix 

(0.005/4  0  \ 

)  (2.6-1) 

0  0.083/ 
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Figure  2 


Figure  2 


Upper-Zone  Tension-Water  Element 


■:5 


?«.  -I - 1 - 

:  5  io  ->* 

Time 


.6-4  Upper-Zone  Tension-Water  Content  for 

Precipitation  Record  BC1 :  True  Values 
(Solid  Line)  and  Predicted  Values  With 
Measurement  C'ov'arianee  R^  (Dashed  Line) 

Upper-Zone  Tension-Water  Element 


: o  z . s 

.6-5  Upper-Zone  Tension-Water  Content  for 

Precipitation  Record  BC1 :  True  Values 
(Solid  Line)  and  Predicted  Values  With 
Measurement  Covariance  R^  (Dashed  Line) 
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Upper-Zone  Free-Water  Element 


T;ME  (aoyS/ 

Figure  2.6-6  Upper-Zone  Free-Water  Content  for 

Precipitation  Record  BC1 :  True  Values 
(Solid  Line)  and  Predicted  Values  With 
Measurement  Covariance  R^  (Dashed  Linel 


Upper-Zone  Free-Water  Element 


L  z  2 . 5 


Figure  2.6-7  Upper-Zone  Free-Water  Content  for 

Precipitation  Record  BC1 :  True  Values 
(Solid  Line)  and  Predicted  Values  With 
Measurement  Covariance  R.;  (Dashed  Line) 
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Lower-Zone  Tension-Water  Element 


Figure  2.6-8  Lower-Zone  Tension-Water  Content  for 

Precipitation  Record  BC1 :  True  Values 
(Solid  Line)  and  Predicted  Values  With 
Measurement  Covariance  R^  (Dashed  Line) 


Lower-Zone  Tension-Water  Element 


'  CC.v  S 

Figure  2.6-9  Lower-Zone  Tension-Water  Content  for 

Precipitation  Record  BtT]  :  True  Values 
(Solid  Line)  and  Predicted  Values  With 
Measurement  Covariance  (Dashed  Line) 
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Figure  2. 


Figure  2 


Lower-Zone  Primary  Aquifer 


"  r\ 


~!ME  ,  3C*S } 

6-10  Lower -Zone  Primary  Aquifer  Content  for 
Precipitation  Record  BC1 :  True  Values 
(Solid  Line)  and  Predicted  Values  With 
Measurement  Covariance  R^  (Dashed  Line) 


Lower-Zone  Primary  Aquifer 

"0-T - 

I  \ 


Q< 


Ml  . 30 v s 

6-11  Lower-Zone  Primary  Aquifer  Content  for 
Precipitation  Record  BCl :  True  Values 
(Solid  Line)  and  Predicted  Values  With 
Measurement  Covariance  R2  (Dashed  Line) 
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Lower-Zone  Supplementary  Aquifer 


Time:  COySj 


Figure  2.6-12  Lower-Zone  Supplementary  Aquifer  Content  for 

Precipitation  Record  BC1:  True  Values  (Solid 
Line)  and  Predicted  Values  with  Measurement 
Covariance  Rj  (Dashed  Line) 


Lower-Zone  Supplementary  Aquifer 


- '/L  ;JCvS 


Figure  2.6-13  Lower-Zone  Supplementary  Aquifer  Content 

for  Precipitation  Record  BC1 :  True  Values 

(Solid  Line)  and  Predicted  Values  With 
Measurement  Covariance  R,;  (Dashed  Line) 
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Figure  2 


Figure  2 
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Time  (says'. 

Excess  of  the  Additional  Impervious  Area 
Content  Over  the  Upper-Zone  Tension-Water 
Content  for  Precipitation  Record  BC1  : 

True  Values  (Solid  Line)  and  Predicted 
Values  with  Measurement  Covariance  R. 
(Dashed  Line) 


TimE  3o.  s 

Excess  of  the  Additional  Impervious  Area 
Content  Over  the  Upper-Zone  I ens i on - Wa t e r 
Content  for  Precipitation  Record  BC 1 :  True 
Values  (Solid  Line)  and  Predicted  Values 
with  Measurement  Covariance  R^  (Dashed  Line) 


2-50 


The  analytic  sciences  corporation 


Six-Hour  Channel- Inflow 


t:m£  (days'. 


Figure  2.6-16  Six-Hour  Accumulated  Channe 1  - 1 n f 1 ow  for 

Precipitation  Record  BC1 :  True  Values 
(Solid  Line)  and  Predicted  Values  With 
Measurement  Covariance  R ^  (Dashed  Line) 


Six-Hour  Channel  -  Inf low 


'  V£  QQ/S- 


Figure  2.6-17  Six-Hour  Accumulated  Channe 1  -  1 n f 1 ow  fo. 

Precipitation  Record  BCl :  True  Values 
(Solid  Line)  and  Predicted  Values  With 
Measurement  Covariance  R^  (Dashed  Line) 
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Discharge  Rate 


Figure  2.6-18  Instantaneous  Discharge  for  Precipitation 

Record  BC1 :  True  Values  (Solid  Line)  and 
Predicted  Values  with  Measurement  Covariance 
R^  (Dashed  Line) 


Discharge  Rate 


ME  OOvSi 


Figure  2  6-19  Instantaneous  Discharge  for  Precipitation 

Record  BC1 :  True  Values  (Solid  Line)  and 
Predicted  Values  With  Measurement  Covariance 
R2  (Dashed  Line) 
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The  root-mean-square  ( rms )  errors  associated  with  R,  are  o,= 

1  3  1 

0.073  mm  for  the  precipitation  measurement  and  02=0.29  m  /sec 

for  the  discharge  rate  measurements.  These  values  correspond 

to  errors  of  round-off  to  the  closest  1/100  of  an  inch  for 

3 

precipitation  and  1  m  /sec  for  discharge.  Odd  numbered  fig¬ 
ures  correspond  to  the  covariance  matrix 


R 


2  ' 


(2.6-2) 


for  which  the  rms  errors  are  five  times  larger  than  those  of 


R 


1 ' 


The  predicted  values  of  instantaneous  channel  dis¬ 
charge  obtained  with  the  use  of  the  covariance  matrices  R^  and 
R2  are  presented  in  Figs.  2.6-18  and  2.6-19,  respect ively . 

The  predicted  discharge  agrees  better  with  the  true  values 
when  the  larger  covariances  are  used.  The  rms  difference  be¬ 
tween  predicted  and  true  discharge  is  7.8  rn  /sec  for  R,  and 

3  J 

7.2  m  /sec  for  R2 .  The  same  type  of  observation  applies  to 

all  states  of  the  SSM  model  (Figs.  2.6-9  through  2.6-17)  ex¬ 
cept  for  the  lower-zone  supplementary  free-water  content  (Figs. 
2.6-12  and  2.6-13)  where  the  fi.  obtained  with  the  use  of  Rj 
is  slightly  better  than  that  obtained  with  R?. 

The  fact  that  larger  measurement  error  covariances 
produce  better  results  can  be  explained  as  follows:  With  a 
small  measurement  error  covariance  the  measurements  are  pre¬ 
sented  to  the  filter  as  being  extremely  precise.  As  a  conse¬ 
quence  small  differences  between  the  predicted  and  measured 
values  are  heavily  weighted  by  the  filter  (i.e.,  t he  Kalman 
gains  are  very  large).  The  filter  responds  quickly  to  dif¬ 
ferences  between  estimated  and  observed  values  and  tends  to 
have  a  fast  oscillatory  response.  This  behavior  is  particu- 
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larly  well  exemplified  in  the  estimates  of  the  upper  and  lower 
zone's  tension-water  content  (Figs.  2.6-4  and  2.6-8).  Note 
that  these  oscillations  happen  to  occur  when  the  upper-zone 
tension-water  element  is  full  and  in  some  instances  (see  Fig. 
2.6-6)  when  surface  runoff  is  being  produced.  When  the  meas¬ 
urement  error  covariance  is  increased,  the  filter  places  more 
trust  in  the  estimates  it  is  computing  and  small  differences 
between  predicted  and  measureed  values  are  given  less  emphasis. 
The  response  of  the  filter  is  somewhat  slower  but  its  overall 
performance  is  improved. 

The  initial  state  covariance  matrix,  P(0),  used  in 
obtaining  the  results  given  in  Figs.  2.6-4  through  2.6-19  was 

taken  to  be  a  diagonal  matrix,  P'  ,  whose  elements  along  the  dia- 

2  2  2 

gonal  were  pjj  =  1.44  mm  ,  p^  =  0.0225  mm  ,  p^  =  2.56  mm  , 

p^  =  1.96  mm^ ,  p^  =  0.0196  nun^,  p^  =  2.56  nim^  and  pj^  =  0 
for  i  >  6.  These  values  correspond  to  an  initial  state  uncer¬ 
tainty  of  1%  of  the  capacity  of  each  of  the  elements  of  the 
SSM  model . 


In  order  to  examine  the  behavior  of  the  filter  for 
different  assumptions  on  the  initial  state  uncertainty,  a  com¬ 
parison  was  made  of  the  results  obtained  for  the  Bird  Creek 
model  with  the  use  of  the  initial  state  covariance  P'  and  those 

obtained  with  a  larger  covariance  matrix,  P" ,  with  diagonal 

2  2  2 
elements  pj'^  =  144  nun  ,  p^  =  2.25  mm  ,  p^  =  256  mm  ,  p^  r 

196  mm^ ,  p^  =  1.96  mm^  ,  p^  =  256  inm^  and  p\’  .  =  0  for  i  >  6, 

corresponding  to  an  rms  uncertainty  ol  10%  of  the  capacity  of 


*The  maximum  value  of  x^,  the  excess  of  the  additional  imper¬ 
vious  area  content  over  the  upper-zone  tension  water  content 
is  x^,  the  capacity  of  the  lower-zone  tension  water  element. 
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the  states  of  the  SSM  model.  Six-hour  measurements  of  pre¬ 
cipitation  and  instantaneous  discharge  were  used  by  the  fil¬ 
ter.  The  measurement  error  covariance  matrix  was  . 

To  prevent  the  state-covariance  from  changing  sig¬ 
nificantly  from  its  initial  values  before  abrupt  changes  in 
the  system  occur,  the  precipitation  record  BC2 ,  shown  in  Fig. 
2.6-20,  was  used  in  the  comparison.  The  results  are  presented 
in  Figs.  2.6-21  through  2.6-36.  Odd  numbered  figures  corre¬ 
spond  to  results  obtained  with  the  initial  state  covariance 
P',  even  numbered  figures  to  those  obtained  with  P" .  As  be¬ 
fore,  true  values  and  one- step-ahead  filter  predictions  are 
indicated  by  solid  and  dashed  lines,  respec t i ve 1 y . 

Six-Hour  Accumulated  Precipitation 


Time  (.days) 


Figure  2.6-20  Precipitation  Record  BC2 
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Figure  2 


F i gu  re 


Upper-Zone  Tension-Water  Element 


.6-21  Upper-Zone  Tension-Water  Content  for 

Precipitation  Record  BC2 :  True  Values 
(Solid  Line)  and  Predicted  Values  With 
P(0)  =  P'  (Dashed  Line) 

Upper-Zone  Tension-Water  Element 


.6-22  Upper-Zone  Tension-Water  Content  for 

Precipitation  Record  BC2 :  True  Values 
(Solid  Line)  and  Predicted  Values  With 
P(0)  =  P"  (Dashed  Line.-) 
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Upper-Zone  Free-Water  Element 


Figure  2.6-23 


TIME  (doys) 

Upper-Zone  Free-Water  Content  for 
Precipitation  Record  BC2 :  True  Values 
(Solid  Line)  and  Predicted  Values  With 
P(0)  =  P'  (Dashed  Line) 


Upper-Zone  Free-Water  Element 
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Figure  2.6-26 


TIME  idovs; 

Upper-Zone  Free-Water  c.ontent  for 
Precipitation  Record  BC2  :  I  rue  Values 

(Solid  Line)  and  Predicted  Values  With 
P(0)  =  P"  (Dashed  Line) 
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Lower-Zone  Tension-Water  Element 


Figure  2.6-25  Lower-Zone  Tension-Water  Content  for 

Precipitation  Record  BC2 :  True  Values 
(Solid  Line)  and  Predicted  Values  with 
P(0)  =  P'  (Dashed  Line) 


Lower-Zone  Tension-Water  Element 


Lower-Zone  lens i on-Wa t e r  Content  for 
Precipitation  Record  BC2 :  True  Values 
(Solid  Line)  and  Predicted  Values  With 
P(0)  =  P”  (Dashed  Line) 


E i gu re  2.6-26 
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Figure  2 


Figure  2 


Lower-Zone  Supplementary  Aquifer 


6-29  Lower-Zone  Supplementary  Free-Water  Content 
for  Precipitation  Record  BC2 :  True  Values 
(Solid  Line)  and  Predicted  Values  With 
P(0)  =  P'  (Dashed  Line) 


Lower-Zone  Supplementary  Aquifer 


for  Precipitation  Record  BC2 :  True  V. tines 

(Solid  Line)  and  Predicted  Values  With 
P(0)  =  P"  (Dashed  Line) 
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Six-Hour  Channe 1  -  I n f 1 ow 


TIME  (doys) 


Figure  2.6-33  Six-Hour  Accumulated  Channel  -  Inf low  for 

Precipitation  Record  BC2 :  True  Values 
(Solid  Line)  and  Predicted  Values  With 
P(0)  =  P*  (Dashed  Line) 


Six-Hour  Channel - In fl ow 


Figure  2.6-34  Six-Hour  Accumulated  Channe J -  1 n f 1 ow  for 

Precipitation  Record  BC2 :  True  Values 
(Solid  Line)  and  Predicted  Values  With 
P(0)  =  P”  (Dashed  Line) 


THE  ANALYTIC  SCIENCES  CORPORATION 


THE  ANALYTIC  SCIENCES  CORPORATION 


Surface  runoff  occurred  on  three  different  occasions: 
on  days  two,  three  and  twenty  (see  Figs.  2.6-23  and  2.6-24). 
Instantaneous  discharge  plots  are  given  in  Figs.  2.6-35  and 

2.6- 36.  For  both  initial  state  covariance  matrices,  the  rms 

differences  between  the  predicted  and  true  discharges  are  very 

3 

close:  7.2  m  /sec  for  the  case  in  which  P'  was  used  and  7.3 

3 

m  /sec  for  the  case  in  which  P"  was  used.  The  one-step-ahead 
prediction  of  the  values  of  the  upper-zone  stales  obtained 
with  the  use  of  P'  and  P"  (Figs.  2.6-21  through  2.6-24)  are 
practically  identical  and  agree  very  well  with  the  true  values. 
However,  for  the  rest  of  the  states  of  the  SSM  model  (Figs. 

2.6- 25  through  2.6-32),  the  agreement  between  predicted  and 
true  values  is  much  better  for  the  estimates  obtained  with  the 
use  of  P',  the  smaller  covariance  matrix.  In  fact ,  the  same 
input  sequences  were  used  in  a  filter  in  which  the  initial  rms 
values  for  the  uncertainty  in  the  values  of  the  states  of  the 
SSM  model  was  set  equal  to  20u/o  of  the  maxima  of  the  state  values. 
This  filter  became  unstable  during  the  fifth  day  of  operation. 

Therefore,  it  is  recommended  that  in  practice,  on¬ 
line  operation  of  the  filter  be  started  when  there  is  little 
uncertainty  as  to  the  correct  values  of  the  states  (e.g.,  dur¬ 
ing  the  dry  season)  or  that  previously  collected  precipitation 
and  discharge  data  be  used  to  carrj'  the  filter  into  steady- 
state.  Another  alternative  is  to  try  to  identify  the  initial 
state.  This  possibility  is  discussed  in  Section  4. 

The  White  River  basin  model  was  used  to  test  the  re¬ 
sponse  of  the  filter  for  a  set  of  parameters  other  than  those 
of  the  Bird  Creek  basin.  The  initial  state  values  in  the  sim¬ 
ulation  and  filtering  algorithms  were  arbitrarily  chosen  as 
Xj(0)  =  40  mm,  X2<0)  =  20  mm,  x^(0)  =  200  mm,  x^(0)  =  40  mm, 

Xr}(0)  =  100  mm  and  x.(0)  =  0  for  i  >  5.  For  the  filter,  the 
rms  uncertainty  in  the  initial  states  was  taken  to  be  1%  of 
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the  element's  capacity  for  the  first  six  states  and  zero  for 
the  remaining  states. 

Measurements  produced  by  the  simulation  were  precipi 
tation  at  six-hour  intervals  and  mean  discharge  at  24-hour 
intervals.  Thus,  k=6,  £=1  and  m=4 .  The  prec i pi ta t i on  record 
produced  by  the  simulation  is  depicted  in  Fig.  2.6-37.  The 
measurement  error  covariance  matrix  used  by  the  filter.  R^, 
was  that  implied  by  the  covariance  matrix  R9  for  mean  daily 
discharge ,  i . e . , 

(0.13 

0 


(2.6-3) 


0.52, 


Six-Hour  Accumulated  Precipitation 


Figure  2.6-37 


Precipitation  Record  for  White  River  Basin 
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The  results  obtained  are  presented  in  Figs.  2.6-37 
through  2.6-45.  The  excellent  fit  between  the  one - s t ep- ahead 
(six  hours)  predicted  mean  daily  discharge  and  the  true  dis¬ 
charge  shown  in  Fig.  2.6-45  is  typical  of  the  results  obtained 
for  both  the  Bird  Creek  and  White  River  basins  when  there  is  no 

surface  runoff.  The  mis  difference  between  the  predicted  and 

3 

true  mean  daily  discharge  is  only  1.4  m  /sec. 

It  was  not  possible  to  obtain  simulation  results  for 

the  White  River  basin  in  which  surface  runoff  appeared.  Note 

that  even  with  the  substantial  amount  of  rainfall  produced 

/V 

during  the  last  four  days  of  simulation,  the  drainage  and 
percolation  rates  from  the  upper-zone  free-water  element  (Fig. 

2.6- 39)  preclude  the  value  of  the  state  from  reaching  its  maxi¬ 
mum  and,  thus  prevent  the  appearance  of  surface  runoff. 

The  same  random  sequence  that  was  used  in  the  simu¬ 
lation  to  produce  the  precipitation  record  of  Fig.  2.6-37  and 
the  mean  daily  discharge  indicated  by  the  solid  line  in  Fig. 

2.6- 45  was  used  in  a  simulation  in  which  the  precipitation 
record  of  Fig.  2.6-37  was  reproduced  but  values  of  instantan¬ 
eous  discharge  instead  of  mean  daily  discharge  were  generated 
(k=6,  £= 1,  m=l).  These  results  served  as  measurements  for  a 
Kalman  filter  in  which  the  initial  covariance  matrix  was  t he 
same  as  in  the  previous  case  and  the  measurement  error  covari¬ 
ance  matrix  was  R^ .  The  results  for  the  instantaneous  dis¬ 
charge  rate  prediction  are  given  in  Fig.  2.6-46.  The  rms  dif¬ 
ference  between  predicted  and  true  discharge  is  3.8  m  /sec . 


*The  expected  monthly  precipitation  is  42  mm.  Total  precipi¬ 
tation  for  the  record  of  Fig.  2.6-37  is  86  nun.  For  the  last 
four  days,  the  amount  of  precipitation  is  53.7  mm. 
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Lower-Zcne  Primar 


M 


THE  ANALYTIC  SCIENCES  CORPORATION 


Lower-Zone  Supplementary  Aquifer 


Figure  2.6-42  White  River  Basin.  Lower-Zone  Suppl emen t a ry 

Free-Water  Content:  True  Values  (Solid  Line) 
and  Predicted  Values  (Dashed  Line) 
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Impervious  Area  Content  Over  t  hi'  Fpper-Z.one 
Tension-Water  Content:  True  Values  (Solid 
Line)  and  Predicted  Values  (Dashed  Line) 


Figure  2.6-43 
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Discharge  Rate 


Figure  2.6-46  White  River  Basin.  Instantaneous  Discharge: 

True  Values  (Solid  Line)  and  Predicted  Values 
(Dashed  Line) 


The  results  obtained  show  that  notwithstanding  the 
strong  nonlinearities  of  the  soil  moisture  model,  the  extended 
Kalman  filtering  technique  described  in  Section  2 . S  performs 
extremely  well  on  the  combined  basin  model.  The  occurrence  of 
the  peaks  of  the  predicted  discharge  and  their  magnitude  agree 
very  well  with  the  simulation  results. 
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3.  STATE-SPACE  MODEL  DEVELOPMENT  FOR 

UNIT  HYDROGRAPHS 


The  objective  of  the  work  performed  under  this  task 
was  to  develop  a  computer  program  which  will  create  state- 
space  models  approximating  the  impulse  response  of  a  unit  hy¬ 
drograph.  This  computer  program  is  intended  to  be  included  in 
the  NWSRFS  forecast  component.  A  users  manual  documents  the 
computer  program  and  its  operation  in  detail  (Ref.  1).  The 
computer  program  has  the  following  capabilities. 

•  Creation  of  a  d i sc  ret e - t ime  state-space 
model  (with  selectable  output  time  step) 
based  on  the  input  of  a  unit  hydrograph 

•  Production  of  an  output  file  of  river 
discharge  predictions  calculated  from 
the  state-space  model  and  an  input  file 
of  channel  inflow  values 

•  Output  of  graphical  summaries  and  tabula¬ 
tions  of  the  state-space  model  impulse 
response,  and  the  squared  magnitude  spec¬ 
trum  and  phase  spectrum  of  the  transfer 

f unc  t i on . 

The  canonical  variate  method  is  used  in  the  approxi¬ 
mation  of  unit  hydrographs  by  reduced-order  models  ol  state- 
space  form.  This  section  describes  the  mathematical  method  by 
which  the  state-space  models  are  created  and  presents  some 
results  based  on  unit  hydrographs  supplied  by  NOAA/NWS  to  TASC 
for  use  in  this  study.  Section  3.1  formulates  the  problem  of 
unit  hydrograph  approximation  as  a  reduced-order  filtering 
problem.  Section  3.2  presents  the  concept  of  cationie.il  variate 
decomposit  ion  of  the  past  and  future  of  a  random  process.  the 
optimality  of  several  canonical  variate  procedures  are  discussed 
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in  Section  3.3.  The  decomposition  of  past  and  future  is  used 
in  Section  3.4  to  compute  explicitly  the  optimal  reduced-order 
state-space  model.  Section  3.5  presents  results  of  Lhe  appli¬ 
cation  of  this  method  to  unit  hydrograpus  supplied  by  the  NWS. 


3.1  FORMULATION  OF  THE  PROBLEM 

Consider  the  problem  in  which  a  unit  hydrograph  h(i) 
is  given  that  specifies  the  response  at  lag  t  to  a  unit  input 
at  time  zero.  It  is  desired  to  find  a  state-space  model,  pref¬ 
erably  of  low  order,  which  is  a  good  approximat ion  in  some 
sense  to  the  given  unit  hydrograph.  This  problem  cannot  be 
separated  from  the  characteristics  of  the  input  process  since 
the  modes  of  h(t)  that  are  excited  and,  hence,  the  output  de¬ 
pend  strongly  upon  the  input  process.  Nominally,  it  will  be 
assumed  that  the  input  process  is  white  noise  which  excites 
all  frequencies  proportionately.  If  the  typical  input  signal 
power  spectrum  is  known  and  different  from  white  noise,  this 
fact  can  be  easily  included  in  the  method  described  below  and 
would  lead  to  an  alternative  approximating  state-space  model. 

It  will  be  shown  in  Section  3.5  ,hat  the  white  noise  assumption 
leads  to  excellent  approximations  of  the  unit  hydrograph  with 
low  order  state-space  models. 

A  schematic  description  of  the  problem  is  shown  in 
Fig.  2.1-1.  The  problem  of  determining  a  state-space  model 
which  does  a  "good"  job  of  predicting  the  output  v(t)  from  the 
input  u(t)  can  be  viewed  as  a  reduced-order  filtering  problem. 
Consider  u ( t )  and  v(t)  as  two  related  random  processes.  (In 
Section  3.4.3  we  will  consider  the  case  of  a  vector  output 

process  v ( t  )  .  )  Given  the  past  of  u ( t  )  for  t  =0,-1 .  it  is 

desired  to  predict  the  future  evolution  of  v(t)  for  t=l,2«...  . 

A  recursive  or  state-space  filter  of  some  specified  state  order 
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R — 4  596* 

v(t)  OUTPUT 


v(t)  ESTIMATE 


Figure  3.1-1  Approximation  of  Unit  Hydrograph 
by  a  Reduced  Order  Filter 

will  be  derived  from  the  unit  hydrograph.  The  approach  anti 
criterion  of  approximation  is  described  in  the  next  two  sec¬ 
tions,  and  the  filter  derived  in  the  following  section. 

3.2  CANONICAL  VARIATE  DECOMPOSITION  OF  PAST  .AND  FUTURE 

The  central  concept  in  the  approach  involves  use  of 
the  canonical  variate  decomposition  of  the  past  of  one  random 
process  and  the  future  of  another  process  (Refs,  9,  10,  and 
11).  This  corresponds  ter  a  particular  coord  i  nat  i  za  t  i  on  of  the 
predictor  space  (Ref.  12)  in  a  way  that  leads  very  natural 'y 
to  the  selection  of  reduced-order  models  and  filters  which  are 
optimal  in  several  senses  discussed  in  Section  3.3. 

Suppose  u(t)  and  y(t),  for  t  an  integer,  are  time 
series  which  are  jointly  stationary  in  the  wide  sense.  Ve  ore 
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primarily  interested  in  prediction  of  the  future  evolution  of 
y(t)  based  upon  the  past  realization  of  u(t)  and  so  consider 
the  two  vectors  of  random  variables 


1  u(t>  \ 
f  U(t-l) 

/v  (  t  +  1  )  \ 

y  ( t+2  )  ' 

Zj ( t)  = 

iu(  t-p+1  )J 

,  z2( t  )  = 

lv(t+p )) 

(3.2-1  ) 

where  p  is  some  specified  number  of  shifts.  We  call  z^(t)  the 
past  of  u(t)  and  z2(t)  the  future  of  v(t).  The  integer  p  is 
generally  somewhat  larger  than  the  maximal  "state"  order  to  be 
considered  but,  in  theory,  may  be  infinite. 


In  the  reduced-order  filtering  problem,  the  relation¬ 
ship  between  the  past  z ^ { t )  and  the  future  z 2 (t)  is  to  be  ap¬ 
proximated  by  a  Markov  process  model  of  specified  state  order 
k.  In  particular,  the  k^-order  state  expressible  as  k  linear 
combinations  of  the  past  z^(t)  that  best  predicts  the  future 
z2(t)  is  to  be  determined.  For  the  immediate  discussion  "best" 
means  best  percent  error  prediction  which  corresponds  to  maximal 
reduction  in  error  in  prediction  of  z2(t)  on  a  percent  error- 
basis  as  given  by  the  canonical  correlation  criterion  discussed 
in  Section  3.3.  The  best  k  linear  combinations  of  /,  j  (  t  )  for- 
predicting  z2(t)  are  those  having  maximal  correlation  with 
z2(t).  Finding  these  best  linear  combinations  of  z ^ ( t  )  hav  ing 
maximal  correlation  with  z2(t)  is  precisely  the  classical  canoni¬ 
cal  correlations  and  variates  problem  of  mathematical  statistics 
(Refs.  13,  14).  The  more  genera]  canonical  prediction  criterion 
is  discussed  in  Section  3.3.2  as  a  simple  modification  of  the 
classical  canonical  correlations  and  variate  problem.  The 
solution  to  the  canonical  variate  problem  is  obtained  by  put¬ 
ting  the  covariance  structure  of  Zj(t)  and  z.;(t)  in  a  canonical 
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form.  Nonsingular  transformations,  J  and  L,  of  the  past 
and  future  z^  to  new  sets  a  and  b  of  canonical  variables  or 
variates 


a  =  Jz2,  b  =  lz2  (3.2-2) 

are  to  be  found  such  that  in  this  new  basis  the  correlations 
between  the  past  a  and  future  b  have  a  canonical  structure 


cov(a,a)  =  1,  cov(b,b)  =  I 


/ 


cov(a,b)  = 


\ 


0 


0  ^ 

(3.2-3) 

0 


with  the  canonical  correlations  r^  >_  .  .  .  >_  r£  >  0  in  descending 
order.  Thus,  the  components  of  the  past  a  are  mutually  uncor¬ 
related  as  are  those  of  the  future  b.  Of  all  linear  combina¬ 
tions  of  Zj  and  ^ .  the  first  component  of  a  has  maximum  cor¬ 
relation  with  the  first  component  of  b.  It  can  be  shown,  for 
any  order  k,  that  the  first  k  components  of  a  (i.e.,  corre¬ 
sponding  linear  combinations  of  the  past  z.j)  lead  to  the  best 
prediction  of  the  future  z.^.  The  canonical  correlations 
r^+ j  ,  .  .  .  , r^  corresponding  to  the  neglected  variables  give  a 
measure  of  the  amount  of  information  lost  in  using  k  rather 
than  i  components.  The  requirements  of  3.2-3  are  equivalent 
to  finding  J  and  L  such  that 
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J 


cov(z^ , Z2 ) L 


T 


0 


\° 


J  COV(Zj 


z:  )J 


=  I 


L  cov(z2  ,  Z2 ) L1  =  I 


(3.2-4) 


(3.2-5) 


This  is  easily  accomplished  using  a  singular  value  decomposi¬ 
tion  (Ref.  15)  which  is  computationally  very  efficient  and  nu¬ 
merically  very  accurate  and  stable.  Dimensions  of  z^(t)  and 
^(t)  as  high  as  several  hundred  can  be  handled  efficiently 
and  accurately  using  these  computational  techniques. 


3.3  OPTIMALITY  OF  CANONICAL  VARIATE  ANALYSIS 

The  canonical  variate  analysis  described  above  provides 
an  optimal  choice  of  a  restricted  number  of  random  variables 
from  one  set  for  prediction  about  a  second  set  of  random  vari¬ 
ables.  As  classically  formulated  in  Ref.  13,  this  involves  a 
canonical  correlation  criterion  which  is  optimal  in  the  sense 
of 


•  Maximizing  correlation  between  the  ob¬ 
served  set  and  the  predicted  set  of 
random  variables 

•  Maximizing  the  mutual  information  be¬ 
tween  the  chosen  variates  used  for  pre¬ 
diction  and  the  variables  predicted. 
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A  new  generalization,  the  canonical  prediction  criterion,  is 
also  discussed  below  which  is  optimal  in  the  sense  of 

•  Minimizing  the  expected  weighted  squared 
prediction  error. 

The  canonical  prediction  criterion  applied  to  unit  hydrographs 
gives  much  lower  order  state-space  models  which  adequately 
approximate  the  unit  hydrographs.  Comparisons  of  state-space 
models  fitted  using  both  the  canonical  correlation  and  pre¬ 
diction  criteria  are  given  in  Section  3.5.1. 

3.3.1  Canonical  Correlation  Criterion 

The  canonical  variate  problem  was  originally  formu¬ 
lated  (Ref.  13,  also  see  Ref.  14)  as  a  sequential  selection 
procedure.  As  discussed  in  Section  3.2,  the  procedure  is  con¬ 
cerned  with  the  optimal  selection  of  k  linear  combinations  of 
a  vector  z^  of  random  variables  for  optimal  prediction  of  a 

related  vector  z?  of  random  variables.  First  a  pair  of  linear 
z  t  t 

combinations  a^  =  and  b^  =  £^z ^  of  the  respective  vectors 

of  random  variables  z^  and  z^  are  determined  wThich  have  maxima] 
correlation.  Next,  a  second  pair  ~  2-2-1  anc*  ^2  =  -2-2  are 
found  which  are  uncorrelated  with  a^  and  b^  and  which  have 
maximal  correlation.  The  procedure  continues  up  to  the  speci¬ 
fied  number  k  of  linear  combinations  of  z^  which  are  permitted. 
Thus  the  canonical  variate  procedure  finds  the  k  mutually  uncor¬ 
related  components  of  z^  which  are  maximally  correlated  with 
Z2 •  Hotelling  (Ref.  13)  defines  an  intuitive  scalar  measure 
of  correlation  between  the  two  vectors  of  random  variables  Zj 
and  £2 ,  called  the  vector  correlation  coefficient,  which  is 
simply  expressed  in  terms  of  the  non-zero  canonical  correla¬ 
tions  r,  ,  .  .  .  ,  r  as 
1  n 
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Q2  =  r2r2. . .rn  (3.3-1) 

It  can  be  shown  that  the  choice  of  k  linear  combinations  of 
which  have  maximal  vector  correlation  coefficient  are  the  first 
k  canonical  variates. 

The  canonical  correlation  method  maximizes  the  mutual 
information.  Shannon  and  Weaver  (Ref.  16)  define  the  informa¬ 
tion  in  one  random  vector  z about  another  random  vector  z2 , 
now  commonly  called  mutual  information,  by 

r  p 1 2  i ’ —2  ^ 

J(-l;-2}  =J P12(51.z2)  log  p  1  (z  ^ )  p2(z2)  d?l  d?2 

(3.3-2) 

where  the  base  of  the  logarithm  is  arbitrary  and  determines 
the  particular  units  of  information,  and  where  Pj9  is  the  joint 
and  p^  and  p2  the  marginal  probability  densities.  Gelfand  and 
Yaglom  (Ref.  17)  showed  that  the  mutual  information  is  simply 
expressed  in  terms  of  the  canonical  correlations  rj,...,rn 
between  the  two  vectors  by 

J(z1;z2)  =  -  2  X)  log  (1  -  r2 )  =  -  2  l°g  w  (3.3-3) 

where  Hotelling  (Ref.  13)  defines  the  vector  alienation  coef- 
f icient 

w  =  (1  -  r2)  . . .  (1  -  r2)  (3.3-4) 

as  a  measure  of  independence  of  z^  and  z2 .  Gelfand  and  Yaglom 
(Ref.  17)  extend  the  definition  of  mutual  information  to  vectors 
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of  infinitely  many  random  variables,  e.g.,  random  processes  in 
both  continuous  time  and  discrete  time.  This  development  also 
provides  the  basis  for  extending  canonical  variates  to  random 
processes  (Ref.  18). 


Now,  if  a  restricted  number  k  of  linear  combinations 
(a^,..,a^)  of  one  random  process  u(t)  are  used  to  predict 
another  random  process  v(t),  then  the  choice  maximizing  the 
mutual  information  is  the  first  k  canonical  variates  and  the 
mutual  information  is  expressed  by  the  first  k  canonical  cor¬ 
relations 


=  “  \  l°g  H  (1  -  r2.) 

1  j  =  l  J 

(3.3-5) 

Thus  the  canonical  correlation  method  provides  an  optimal  pro¬ 
cedure  in  terms  of  mutual  information  for  choosing  a  finite 
state  representation  of  one  random  process  for  prediction  of 
another . 


max 


J(a1 


.  ak ;v( t ) ) 


3.3.2  Canonical  Pred i c t i on  Criterion 

A  more  general  criterion  of  prediction  error  is  the 
expected  weighted  square  prediction  error 

h(^2  “  | 2 )  =  E(z2  -  I2 )T  6  *(z7  (3.3-b) 

where  0  is  an  arbitrary  positive  definite  symmetric  weighting 
matrix,  is  the  minimum  variance  estimate  of  z9  based  upon 
the  k  selected  linear  combinations  of  Zj  ,  and  E  is  expectation. 
The  optimal  choice  of  k  linear  combinations  of  z^  that  mini¬ 
mizes  h  corresponds  to  the  first  k  canonical  variates  for  a 
canonical  correlation  analysis  with  the  "correlation"  structure 
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(3.3-7) 


although  this  is  not,  in  general,  a  covariance  matrix.  That 
is,  the  true  covariance  Z  22  -2  *s  rePlacec^  by  the  weighting 

matrix  0.  The  minimum  prediction  error  is  simply  expressed  by 


lin  h(z2  -  Z2 )  =  tr 


-  (rl  + 


+  rk> 


(3.3-8) 


in  terms  of  the  canonical  predictors  r^  which  play  the  same 
role  as  the  canonical  correlations. 


Thus  is  can  be  seen  that  when  0  is  set  equal  to 
the  canonical  correlation  problem  can  be  viewed  as  weighting 
the  squared  error  by  the  inverse  covariance  so  that  the  per¬ 
cent  error  or  error  relative  to  the  variance  of  each  variable, 
is  the  criterion  considered.  The  criterion  given  by  Eq .  3.3-7 
is  more  general  than  that  usually  considered  in  the  canonical 
variate  method  and  permits  arbitrary  quadratic  weighting  of 
the  prediction  errors. 

Such  weighting  is  particularly  useful  in  reduced-order 
state-space  modeling  of  unit  hydrographs  and  permits  weighting 
the  more  important  variables  to  be  estimated.  In  particular, 
it  is  found  that  the  criterion  of  the  expected  sum  squared 
prediction  error 


h(z.2  -  z^)  =  E(z2  -  I2'  (-2  ’  - 2 ^ 


(3.3-9) 


of  the  predicted  future  gives  markedly  improved  state-space 
model  approximat ions  to  some  unit  hydrographs.  This  corre¬ 
sponds  to  setting  0  equal  to  I,  the  identity  matrix. 


J 
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3.4  OPTIMAL  REDUCED-ORDER  MODELING 

The  canonical  variates  and  correlations  analysis  is 
used  in  this  section  to  derive  optimal  reduced-order  filters. 

In  the  reduced-order  filter  problem,  a  process  u(t)  and  a  re¬ 
lated  process  v(t)  are  given,  and  predictions  of  v(t)  based 
upon  past  observations  u(t),  u(t-l),...  using  a  reduced-order 
filter  are  required.  As  in  the  canonical  variate  analyis  in 
Section  3.2,  let  z^  ( t )  be  the  past  of  u(t)  and  z^t)  be  the 
future  of  v(t)  (e.g.,  Eq .  3.2-1). 

3.4.1  State  Vector  Determination 

For  a  given  order  k  for  a  reduced-order  filter,  find¬ 
ing  a  best  k-element  state  is  equivalent  to  finding  the  k  linear 

T  T  T 

combinations  of  the  past  z,  =  (u  (0),u  (-1),...)  which  have 

A  rp  rp  rp 

the  best  ability  to  predict  the  future  z^  =  (v  ( 1 ) ,v* (2 ),...) . 

If  the  best  predictive  ability  means  to  minimize  percent  error 
in  predicting  all  components  of  ,  then  this  problem  is  pre¬ 
cisely  the  canonical  variates  and  correlations  problem  of  Section 
3.3.1.  Or  more  generally  we  can  use  the  canonical  prediction 
criterion  of  Section  3.3.2.  To  solve  this  problem,  nonsingular 
transformations  of  the  past  z^  and  future  ^  are  determined 

a  =  JZj ,  b  =  LZ2  (3.4-1) 

such  that  in  this  new  basis  the  past  a  and  future  b  have  a  ca¬ 
nonical  structure  as  in  (3.2-3). 

Once  the  canonical  variate  problem  is  solved,  an  op¬ 
timal  reduced-order  filter  of  order  k  can  be  determined  for 
each  k  <$.  with  the  mi  n  imal -order  realization  given  for  k  =  i!  .  If 
k  linear  combinations  of  the  past  z^  are  t o  be  used  to  predict 
the  future  ^ ,  then  the  optimal  choice  is  the  first  k  canonical 
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variates.  These  k  linear  combinations  specify  the  optimal 
state  vector  of  the  k^-order  filter.  Specifically,  the  state 
vector  is 

x(t)  =  Jkz1(t)  where  Jk  =  (Ik,0)J  (3.4-2) 

with  lk  the  k*k  identity  matrix. 

3.4.2  State-Space  Realization 

The  remainder  of  the  problem  is  to  obtain  a  state- 
space  realization  of  the  optimal  reduced-order  filter.  In 
particular,  expressions  for  the  calculation  of  the  state-space 
matrices  in  terms  of  the  canonical  variate  analysis  are  needed. 
The  desired  state-space  form  of  the  filter  is 

x(t)  =  4>x  ( t  - 1  )  +  Gu(t )  (3.4-3) 

v ( t )  =  Hx(t-l)  (3.4-4) 

with  the  output  v(t)  the  optimal  reduced-order  filter  estimate 
of  v(t).  The  output  v(t)  is  related  to  the  state  x(t-l)  at 
time  t-1  to  insure  a  lag  between  the  input  and  output,  i.e., 
so  that  an  input  does  not  produce  an  instantaneous  output. 

This  is  a  conventional  discrete  time  system  description  -  see 
Section  3.4.4  for  the  modi fi cat  ion  involved  when  the  input 
u(t)  is  a  continuous  accumulation  of  channel  inflow  over  the 
time  interval  (t-l,t).  By  calculating  the  covariance  of  x(t) 


with  x ( t - 1 )  and 

u(t) 

using 

Eq .  3.4-3,  the  matrix  <t> 

and  G  must 

satisfy 

2  ( x  ^  ,  x  ^  ^  j  ) 

=  4> 

I(^t-1 

,xt_j  )  +  G 2 ( u t  ,xt  _j  ) 

(3.4-5) 

2  ( x  ,  u  ) 

=  4> 

i(*t-l 

,ut )  +  G  2 (ut ,ut  ) 

(3.4-6) 

3-12 


THE  ANALYTIC  SCIENCES  CORPORATION 


FILTER 


Figure  3.4-1 


Partial  Innovations  Represent  at i on  of 
Optimal  Reduced  Order  Filter 


where  I(  ,  )  denotes  the  covariance  matrix  between  two  vectors 

and  the  shorthand  xt  -  x(t)  is  used. 

These  equations  are  easily  solved  for  $  and  G  explicit 
ly  in  terms  of  the  various  covariance  matrices.  The  filter 
(3.4-3)  and  (3.4-4)  can  be  put  into  feedback  form 


x ( t )  =  A  x(t-l)  +  K(u(t)  -  C  x ( t - 1 ) ) 


(3.4-7) 


where  the  matrices  A,  K.  and  C  are 


A  =  Mxt.xt_,)  1 


(3.4-8) 


3-13 


THE  ANALYTIC  SCIENCES  CORPORATION 


K  =  U(xt,u£)  -  ^(xt,xt_j)  l"1  (x,._1  ,xt_1 )  2(xt_j,ut)] 
U(ut.ut)  -  1  (ut  .xt_i )  Z-1  (xt_j  ,xt_j  )  Z(xt_j  ,ut )  )_1 

(3.9-9) 

C  =  Z(ut,xtl)  Z_1 (xt _j  ,xt_j  )  (3.9-10) 

The  state-space  form  (3.9-3)  through  (3.9-9)  has  matrices  <t> 
and  G  given  by 

4>=A-KC,G  =  K  (3.9-11) 

These  matrices  have  a  simple  interpretation  as  regression  coef¬ 
ficients.  The  matrix  A  is  the  regression  of  x(t)  on  x(t-l). 
i.e.,  the  best  prediction  of  x(t)  given  x(t-l).  The  term  C 
x(t-l)  is  the  best  prediction  of  u(t)  based  upon  x(t-l).  The 
feedback  gain  matrix  K  is  the  regression  of  the  conditional 
random  variables  x(t)  given  x(t-l)  on  the'  conditional  random 
variable  u(t)  given  x(t-l).  This  gives  the  best  prediction 
of  x(t)  based  upon  u(t),  hav'ing  already  accounted  for  the  de¬ 
pendence  upon  x(t-l).  The  quantity  u(t)  -  c  x(t-l)  can  be 
called  the  partial  innovations  of  the  process  u(t)  with  re¬ 
spect  to  the  past  of  the  process  x(t-l)  (i.e.,  the  new  infor¬ 
mation  in  u(t)  uncorrelated  with  the  past  reduced-order  states 
x(t-l),  x(t-2 ),...).  The  remaining  quantity  to  be  specified 
is  the  prediction  gain  matrix 

H  =  i<vt,xt.1)  Z'^Xj.j  ,xt_2)  (3.9-12) 

for  predicting  y(t)  based  upon  x(t-l).  The  partial  innovations 
realization  of  the  optimal  reduced-order  filter  is  illustrated 
in  Fig.  3.9-1  and  has  the  same  structure  as  the  Kalman  filter 
in  terms  of  optimal  extrapolation,  prediction  of  the  measure¬ 
ment,  calculation  of  the  partial  innovation  and  update  (Rel. 
19).  The  matrices  can  be  explicitly  computed  using  Eq .  3.9-2, 
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the  covariance  structure  between  ( t  )  and  ^(t),  and  the  ca¬ 
nonical  form,  Eq .  3.2-3. 

3. A. 3  Multirate  Unit  Hydrographs 

There  are  a  number  of  situations  in  which  the  channel 
inflow  is  the  result  of  accumulations  of  water  over  a  time 
interval  greater  than  the  time  between  discharge  measurements. 
For  example  accumulated  channel  inflow  may  be  calculated  every 
12  hours  corresponding  to  the  precipitation  measurements  and 
channel  outflow  is  measured  every  6  hours.  To  accommodate 
this  situation,  a  multirate  state  space  model  is  needed  which 
has  inputs  every  12  hours  and  outputs  every  6  hours.  By  prop¬ 
erly  setting  up  the  unit  hydrograph  as  a  single  input  multi¬ 
output  unit  hydrograph,  the  above  algorithms  apply  exactly  as 
they  are. 


We  make  the  following  conventions.  Let  u(t)  be  the 
input  channel  inflow,  t=l,2,...,  with  the  input  sample  rate 
normalized  to  1  time  unit.  Let  y(t),  t  an  integer,  be  the 
vector  of  output  channel  outflow  occurring  at  times  (s-li  <  s 
<  (s|  where  ( s ]  denotes  the  least  integer  greater  than  or  equal 

to  s.  Figure  3.4-2  gives  an  illustration  of  the  input,  the 
output  grouping  and  labeling,  the  unit  hydrograph  and  the  pulse 
input  used  in  its  definition. 


If  lit  is  the  output  sample  rate  then 


v(t  ) 


/output  ( t  -  1  +  A  t  ) 


\  out  put  ( t  -  1  +  rM  )( 


(3.4-13) 


where  r  =  1/^t.  The  vector  unit  response  function  g(T), 
t=l . q,  is  defined  in  terms  of  the  unit  hydrograph  h(t)  by 
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Figure  3.4-2  Multirate  Unit  Hydro graph 


(h ( t - 1  +  At) 

h(t-l  +  rAt)^ 


(3.4-14) 


Once  the  conventions  of  v(t)  as  a  vector  output  and 
g(t)  as  a  vector  unit  hydrograph  are  made,  the  canonical  vari¬ 
ate  analysis  and  reduced  order  state  space  models  are  derived 
as  in  Section  3.4.2.  Note  that  the  state  space  model  operates 
at  the  input  sample  rate  so  that  the  state  vector  only  changes 
when  there  is  a  new  input. 


The  transfer  function  is  a  useful  tool  in  studying  a 
multirate  discrete  time  system.  It  can  be  used  to  describe 
the  response  of  the  system  to  a  sampled  sine  wave  at  the  input  . 
In  Section  3.5.2,  an  example  will  show  that  situations  arise 
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where  the  response  to  a  sampled  si..e  wave  input  at  a  given 
frequency  may  include  a  substantial  sampled  sine  wave  compon¬ 
ent  at  much  higher  frequency  than  the  input.  While  this  is  a 
difficulty  in  the  present  forecasting  system  due  to  multirate 
unit  hydrographs,  it  is  completely  avoided  when  using  the  con¬ 
tinuous  catchment  model  as  in  Section  2. 

Consider  the  multirate  system  where  u(t),  t  = 

. . . -2r , r ,o  ,  r  ,  .  .  .  ,  u^  =  u(kr)  is  the  input  and  v(t),  t  = 

-2  ,  - 1 , 0 , 1 , 2  ,  .  .  .  is  the  output  with  r  an  integer.  Let  k  and  £ 
be  integers  so  that  the  output  at  time  t  =  (k-l)r+£,  0  <  £  <  r, 
is  related  to  the  input  by  the  finite  unit  hydrograph  function 
gjgCt)  as 

q 

v(kr+£-r)  =  ^  g^ ( t )  u(kr-tr)  (3.4-15) 

To  find  the  transfer  function,  we  compute  the  Fourier 
transform  V(uj)  of  v(t)  for  -n  <  u>  <  n 

V(m>  =  £  v(t)  e" *wt  =  £  £  v ( kr+f - r )  e‘ im ( kr+t ' r > 

t  =  -oo  £  =  1  k= 


r 

=  L 

£  =  1 


g/T) 


-  itu  rk 
e 


=  L  elw(r"£)  G.(mr)  £  U(u»r)  (3. 4-lb) 

£  =  1  *  r 

where  G^(\)  and  U(A)  are  the  Fourier  transforms  of  g^(k)  and 
u^  respectively. 
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The  transform  U(A),  -  n  <  A  £  n  ,  of  the  input  process 
=  u(kr)  is  periodic  with  period  2 n  and  in  (3.4-lb)  repeats 
up  to  frequency  A  =  ±nr.  This  is  the  Fourier  extension  or 
Nyquist  interpolation  of  the  input  to  a  higher  sample  rate. 

The  transfer  function 


C(u»)  =  Z  eiu,(r‘£)  G  (wr ) 
£  =  1 


(3.4-17) 


is  then  applied  to  this  interpolated  input  function  to  obtain 
the  output  function.  If  the  input  is  a  sine  wave  =  sin(Ak+o) 
for  0  <  A  <  n  ,  then  the  output  time  function  is 


v(  t ) 


E 

uteQ 


C(w) 


1 

r 


sin(  u»t  +  arg  G(ui)  +  0  ) 


(3.4-18) 


where  fi  is  the  periodic  repetition  of  A  in  terms  of  the  output 
sample  rate 


Q 


{oj:w  =  -  ± 
r 


2nm 

r 


uj  <  n  ,  m  integer} 


(3.4-19) 


3.4.4  Continuous  Accumulation  of  Channel  Inflow 

There  is  a  minor  inconsistency  between  the  usual  way 
of  specifying  a  discrete  time  and  a  continuous  time  system  in 
terms  of  the  unit  pulse  response  and  that  of  implementing  a 
continuous/discrete  time  filter.  This  requires  a  minor  modi¬ 
fication  of  the  reduced  order  filter  given  by  Eqs.  3.4-3  and 
3.4-4. 


Figure  3.4-3(a)  shows  a  constant  pulse  input  and  the 
resulting  continuous  unit  pulse  response  in  (b).  A  discrete 
response  is  shown  in  (d)  which  is  conceptually  considered  as 
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R-59559 


(a)  CONTINUOUS  PULSE  INPUT 


<b)  CONTINUOUS  PULSE  RESPONSE 


(c)  DISCRETE  PULSE  INPUT 


(d)  DISCRETE  PULSE  RESPONSE 


Figure  3.4-3  Relationship  Between  Continuous  and  Discrete 
Time  Unit  Pulse  Response 


the  response  to  an  input  (c)  at  time  t  =  0,  i.e.,  the  pulse 
has  been  accumulated  and  applied  at  the  time  corresponding  to 
the  pulse  beginning.  The  discrete  time  system  is  constrained 
so  that  there  is  no  instantaneous  output  and  hence  no  effect 
occurs  until  the  end  of  the  pulse. 


B 
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The  inconsistency  occurs  when  a  usual  continuous/dis¬ 
crete  time  filter  is  implemented  as  in  Chapter  2.  It  is  natural 
in  such  recursive  schemes  to  associate  the  end  time  of  the 
pulse  with  the  accumulated  channel  inflow  rather  than  the  start 
time  of  the  pulse.  If  this  convention  is  adopted  then  the 
state  space  representation  (Eqs.  3.4-3  and  3.4-4)  is  modified 
to 


x  ( t  )  =  4>  x(t-l)  +  Gu  ( t  )  (3.4-20) 

y ( t )  =  Hx ( t )  (3.4-21  ) 

where  the  quantity  Hx(t-l)  of  Eq .  3.4-4  has  been  shifted  one 
time  step.  This  applies  to  both  the  cases  of  single  output 
and  multirate  output  unit  hydrographs.  Note  the  convention 
that  time  is  normalized  so  that  the  time  between  inputs  is 
unity. 


Associating  the  input  time  with  the  end  time  of  the 
pulse  also  introduces  a  phase  shift  in  the  transfer  function 
(3.4-17)  of  elu,r.  Consequently  the  transfer  function  is  modi¬ 
fied  from  (3.4-17)  to 


G(u) )  =  £  e  G.(uir) 

£-\  * 


(3.4-22) 


3.5  RESULTS 

The  reduced-order  state-space  modeling  described  above 
has  been  applied  to  unit  hydrographs  for  a  number  of  basins 
supplied  by  NWS.  The  character  of  the  reduced-order  models  is 
illustrated  in  Section  3.5.1  and  described  in  more  detail  in 
Refs.  2  and  4.  The  problem  of  spurious  high  frequency  behavior 
inherent  in  some  mult  irate  unit  hydrographs  is  discussed  in 
Sect  ion  3.5.2. 
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3.5.1  Reduced-Order  State-Space  Models 

The  differences  in  reduced-order  models  obtained  from 
the  canonical  correlation  and  canonical  prediction  criteria 
depend  very  strongly  upon  the  spectral  shape  of  the  hydrograph 
transfer  function.  A  striking  comparison  in  fit  using  the  two 
criteria  was  obtained  for  the  Bird  Creek  basin.  The  six-hour 
unit  hydrographs  based  upon  the  input  hydrographs  and  the  ca¬ 
nonical  correlation  procedure  are  shown  in  Figs.  3.5-1  and 
3.5-2  for  4-  and  8-state  models  respect ively .  The  respectively 
squared  magnitude  transfer  functions  are  shown  in  Figs.  3.5-3 
and  3.5-4.  Fits  obtained  using  the  canonical  prediction  proce¬ 
dure  are  illustrated  by  the  unit  hydrographs  in  Figs.  3.5-5 
and  3.5-6  and  by  the  squared  magnitude  transfer  functions  in 
Figs.  3.5-7  and  3.5-8.  Note  that  even  the  8-state  unit  hydro¬ 
graph  from  the  canonical  correlation  procedure  has  a  signifi¬ 
cant  nonzero  tail  whereas  the  4-state  unit  hydrograph  from  the 
canonical  prediction  procedure  produces  an  excellent  fit. 
Figures  3.5-3  and  3.5-4  clearly  illustrate  the  tendency  of  the 
canonical  correlation  procedure  to  fit  all  frequencies  with 
nearly  equal  percent  error,  whereas  from  Figs.  3.5-7  and  3.5-8 
it  is  seen  that  in  the  canonical  prediction  procedure  frequency 
bands  of  highest  energy  are  emphasized.  Thus  for  a  hydrograph 
with  a  large  spectral  peak  and  complicated  spectral  shape, 
i.e.,  requiring  a  high  order  rational  function  for  a  good  ap¬ 
proximation,  the  canonical  prediction  criterion  can  be  expected 
to  excel . 

3 .5.2  Multi  rate  Unit  Hydrographs 

The  reduced-order  modeling  of  multirate  unit  hydro- 
graphs  produces  results  very  similar  to  those  described  in 
Section  3.5.1  for  unit  hydrographs  with  the  same  input  and 
output  rates.  The  major  differences  encountered  for  mult  irate 
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Figure  3.5 


Figure  3.5 


4  STATES 


TIME  (days) 

-1  6-Hour  Unit  Hydrographs  for  Bird  Creek, 
Original  Hydrograph  (Solid  Line)  and 
Fourth-Order  State-Space  Approximation 
Using  the  Canonical  Correlation  Procedure 
(Dashed  Line) 
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TIME  (days) 

-2  6-Hour  Unit  Hydrographs  for  Bird  Creek, 
Original  Hvdrograph  (Solid  Line)  and 
Eight-Order  State-Space  Approximation 
Using  the  Canonical  Correlation  Procedure 
(Dashed  Line) 
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Figure  3.5 


Figure  3.5 


R4I4(» 


-3  Squared  Magnitude  Transfer  Function  of  6-Hour 

Unit  Hydrographs  for  Bird  Creek,  Original  Hydro 
graph  (Solid  Line)  and  Four-Order  State-Space 
Approximation  Using  the  Canonical  Correlation 
Procedure  (Dashed  Line)  M4*4oS 
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-4  Squared  Magnitude  Transfer  Function  of  6-Hour 

Unit  Hydrographs  for  Bird  Creek,  Original  Hydro 
graph  (Solid  Line)  and  Eight -Order  State-Spue 
Approximation  Using  the  Canonical  Correlation 
Procedure  (Dashed  Line) 
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TIME  (days) 

Figure  3.5-5  6-Hour  Unit  Hydrographs  for  Bird  Creek, 
Original  Hydrograph  (Solid  Line)  and 
Fourth-Order  State-Space  Approximation 
Using  the  Canonical  Prediction  Procedure 
(Dashed  Line) 

8  STATES 


6-Hour  Unit  Hydrographs  for  Bird  Creek. 
Original  Hydrograph  (Solid  Line)  and 
Eight -Order  State-Space  Approximation 
Using  the  Canonical  Prediction  Procedure 
( Dashed  Line  ) 


Figure  3.5-6 


THE  ANALYTIC  SCIENCES  CORPORATION 


unit  hydrographs  is  inherent  in  their  multirate  nature  and 
involves  the  production  of  components  at  the  output  which  are 
of  substantially  higher  frequency  than  the  input. 

For  example,  consider  the  unit  hydrograph  with  input 
every  6  hrs  and  output  every  3  hrs  shown  in  Fig.  3.5-9  with 
transfer  function  shown  in  Fig.  3.5-10.  As  expressed  in  Eq . 
3.4-18,  an  excitation  by  a  sampled  sine  wave  at  frequency  f 
cyc/day  for  0  <  f  <  2  produces  at  the  output  the  sum  of  two 
sine  waves 

Aj  sin(2nf  +  0 ^ )  +  A2  sin(2n(4-f)  +  0 ^ 

2  2 

where  A^  and  A ^  come  from  Fig.  3.5-10  at  frequencies  f^  and 
4-f^  cyc/day  respectively.  If  a  0.5  cyc/day  diurnal  component 
is  exciting  the  unit  hydrograph,  then  the  output  is  the  sum  ol 
a  0.5  cyc/day  and  a  3.5  cyc/day  with  amplitude  about  one-eighth 
the  0.5  cyc/day  component.  The  input  and  output  are  illustrated 
in  Figs.  3.5-11  and  3.5-12  where  the  boxes  denote  the  sampling 
times. 


The  problem  of  spurious  high  frequency  components  in 
the  output  of  some  multirate  unit  hydrographs  occurs  in  the 
present  NOAA/NWS  forecast  system  because  the  catchment  model 
operates  at  the  same  rate  as  the  precipitation  measurement. 
This  difficulty  is  completely  avoided  by  implementing  the  con¬ 
tinuous  catchment  model  as  discussed  in  Section  2  so  that  a 
channel  inflow  is  produced  at  the  same  rate  as  the  channel 
d i scha  rge . 
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T I  ML  (.joys) 

Figure  3.5-9  Multirate  Unit  Hydrograph  With  6-Hour  Inputs 
and  3-Hour  Outputs 


0.0  0.5  1.0  1.5  2.0  2.5  3.0  35  4  0 


FREQUENCY  (cyc/doy) 


Figure  3.5-10 


Squared  Magnitude  Transfer  Function  of 
Mult  irate  Unit  Hydrograph  Showing  Output 
Components  in  Response  to  a  0.5  cyc/day  Input 
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Figure  3.5-11  Input  Sine  Wave  and  Sampled  Values 


Figure  3.5-12  Output  Function  With  0.5  cyc/day  and 

3.5  cyc/day  Components 
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•  Computation  of  the  ML  estimates  is  an 

optimization  problem  with  many  different 
ways  available  to  take  advantage  of  sparse 
model  structures. 

The  major  computational  problem  is  to  obtain  maximum 
likelihood  estimates  of  the  model  parameters  using  numerical 
optimization.  The  early  development  of  system  identification 
techniques  was  plagued  by  ill-conditioned  optimization  prob¬ 
lems  caused  by  "nonident i f iable"  parameters.  Models  often 
include  parameters  which  have  little  or  no  effect  upon  the 
measurements.  Although  the  statistical  problem  of  estimating 
such  parameters  is  still  well  defined  and  meaningful,  severe 
numerical  problems  can  arise  in  the  optimization  procedure. 

As  discussed  below,  these  problems  can  be  largely  avoided  by 
inspecting  the  Fisher  information  matrix  which  defines  the 
parameters  or  combinations  of  parameters  that  are  not  identi¬ 
fiable.  The  optimization  is  then  constrained  to  those  param¬ 
eters  which  are  identifiable. 

In  Section  A. 1.1,  a  general  description  of  the  dy¬ 
namical  model  of  a  river  basin  for  the  purpose  of  parameter 
identification  is  given.  A  more  detailed  description  is  given 
in  Section  2  and  previous  progress  reports  (Refs.  2,  A,  5). 
Section  A. 1.2  is  devoted  to  evaluation  of  the  likelihood  func¬ 
tion  using  a  Kalman  filter  and  to  related  computations  used  in 
optimization.  A  detailed  description  of  the  Kalman  filter  for 
the  hydrologic  model  was  given  in  Chapter  2  of  this  report. 

The  parameter  sensitivities  of  the  state  estimate  and  its  error 
covariance  which  are  required  in  the  optimization  are  detailed 
in  Section  A. 1.3.  In  Section  A. 2,  optimization  considerations 
are  discussed  including  a  new  i dent i f iabi 1 i t y  theory  which 
avoids  i 1 1 -condi t ioning  due  to  nonident if iabi 1 ity  of  parameters 
Details  of  a  quadratic  optimization  algorithm  and  statistical 
convergence  criterion  are  described.  Finally  reparamet er iza- 
tions  which  will  accelerate  convergence  are  considered.  In 
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Section  A. 3,  results  from  determining  parameter  identi f iabi 1  i  ty 
and  demonstrating  the  algorithms  on  simulated  precipitation  and 
channel  discharge  data  are  presented. 

A . 1  LIKELIHOOD  FUNCTION  FOR  DYNAMICAL  MODELS 

To  apply  the  maximum  likelihood  method  to  estimate 
the  parameters  of  a  dynamical  hydrological  model  from  observa¬ 
tion  data  requires  a  procedure  for  evaluation  of  the  likelihood 
function  for  such  models.  First,  the  form  of  the  parametric 
dynamical  model  is  described,  and  then  compu t a t i ona 1  procedures 
for  evaluation  of  the  likelihood  function  and  related  quanti¬ 
ties  needed  in  the  optimization  procedure  are  discussed. 

A . 1 . 1  Dynamical  Parametric  Models 

The  dynamical  model  of  a  basin  described  in  Chapter  2 
can  be  regarded  as  parametric  in  the  parameters  0  as  listed  in 
Table  2.2-1.  This  model  is  described  by  the  dynamic  relation¬ 
ships 


x(t)=f(x,  u,  w,  0,  t)  (A. 1-1) 

with  initial  condition  x(t  )  =  x  at  time  t  where  x(l)  is  the 

-  o  — o  o 

state  vector,  u(t)  is  a  known  deterministic  input,  w(t)  is 
white  noise  with  covariance  matrix  Q(t)  and  t  is  time.  The 

initial  condition  x  can  be  considered  as  a  fixed  but  unknown 

-o 

constant  and  included  in  the  vector  of  parameters  to  be  esti¬ 
mated.  The  measurement  model  is  of  the  form 

z  ( t  .  )  =  h  (  x  ,  u  ,  y ,  0  ,  t  )  (  A  .  1  -  2  ) 
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where  y(t)  is  white  measurement  noise  with  covariance  matrix 
R(t).  The  noise  covariance  matrices  Q(t,0)  and  R(t,0)  can 
also  be  considered  as  functions  of  unknown  parameters  to  be 
es  t ima  ted . 

These  equations  and  the  correlation  structure  of  the 

system  noise  terms  imply  a  correlation  structure  among  the  ob- 
T  T  T 

servations  z  =  (Zj,...,z^).  To  determine  maximum  likelihood 
estimates  of  the  parameters  0  given  the  observations  z  requires 
a  procedure  for  evaluating  and  maximizing  the  likelihood  function. 

4.1.2  Likelihood  Function  Calculations 

The  optimization  method  outlined  in  section  4.2  re¬ 
quires,  for  a  specified  value  of  the  parameters  0,  calculation 
of  the 

•  Log  likelihood  function  In  p(z,0) 

•  Gradients  8ln  p(z,0)/86^,  i=l . k 

•  Fisher  information  matrix  with  elements 

-E[82  In  p(z,0  )/80.80j  |  ,  i.,j.  =  l . f 

where  p(z,0)  is  the  probability  density  of  the  measurements  z. 

These  computations  are  performed  by  the  modular  and 
TM 

general  PARAIDE  (parameter  ident i f ica t ion  )  computer  software 
for  maximum  likelihood  parameter  identification.  The  calcula¬ 
tions  are  structured  as  follows  and  described  in  more  detail 
in  following  sections: 

•  Linearize  the  differential  and  measure¬ 
ment  equations  (4.1-1)  and  4.1-2)  u~ing 

PARAIDE  is  a  trademark  of  The  Analytic  Sciences  Corporation. 


THE  ANALYTIC  SCIENCES  CORPORATION 


symbolic  differentiation  to  obtain  the 
state  space  matrices  F,  G,  H  specifying 
the  linear  state  space  equations  about 

each  filtered  state  estimate  x.  ... 

-1/1 

•  Implement  a  Kalman  filter  using  the  line¬ 
arized  state  equation  model  to  propagate 
the  state  estimate  and  its  error  covari¬ 
ance  matrix  P  .  , .  . 

1/1 

•  Calculate  the  innovations  sequence  from 
the  Kalman  filter  and  the  resulting  error 
covariance  matrices.  Evaluate  the  log 
likelihood  function. 

•  Calculate  the  sensitivity  functions 
3F/ae,  3G/ae,  aH/ao,  ag/ao,  aR/ae  by 
further  symbolic  differentiation  of  the 
differential  and  measurement  equations 
(A. 1-1)  and  (A. 1-2). 

•  Implement  a  Kalman-type  algorithm  to 

propagate  3x.  ,./BQ  .  and  3P.  ,./BQ  .  for 
H  H  6  -x/i'  j  i/i'th  J 

each  0j  and  to  evaluate  the  j  gradi¬ 
ent  component  of  the  log  likelihood 
function  3ln  p(z,0)/30j. 

•  Evaluate  the  Fisher  information  matrix 
using  the  gradient  calculations. 


A  method  is  described  for  evaluating  the  joint  like- 

T  T  1' 

lihood  function  p(z,0)  of  the  observations  z  =  (z^ . z^. ) 

for  the  model  described  in  Eqs.  A. 1-1  and  A. 1-2  with  a  parti¬ 
cular  assumed  value  0  for  the  parameters.  The  approach  is  to 
take  advantage  of  the  independence  properties  of  the  innova¬ 
tions  of  a  Kalman  filter. 


Consider  a  value  of  0  for  which  evaluation  of  the 
likelihood  function  is  desired.  If  the  observations  z  were 
produced  by  the  model  Eqs.  A. 1-1  and  A. 1-2  with  true  param¬ 
eters  equal  to  0  ,  then  the  Kalman  filter  based  upon  this 
model  would  produce  an  innovations  sequence  v ^  ‘  ‘  '  ’-N 
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V  . 
—  l 


2  . 
-1 


H  .  x  .  , .  . 
i-i/i -1 


(A. 1-3) 


such  that  is  independent  of  y  ^ 
of  v  . 

-l 


for  i  ^  j  with  the  covariance 


B  . 

l 


H  .  P  .  . .  ,  HT  +  R  . 

l  l/i-l  i  l 


(A.  1-4) 


For  a  linear  model  the  innovations  have  a  Gaussian  distribu¬ 
tion,  and  for  mild  nonlinearities  the  model  is  approximately 
Gaussian.  In  any  case,  the  quadratic  term  in  Eq .  A. 1-5  gives 
a  weighted  squared  measure  of  the  innovations  vo  .  Assuming  a 
Gaussian  distribution  for  the  innovations  leads  to  the  expres¬ 
sion  for  the  logarithm  of  the  probability  density 


1  rP  1 

In  p(z.0  )  =  "  T  (In  B.  +  v  1  B.  v'.)  +  constant 

1-1  (A. 1-5) 

where  the  constant  depends  only  on  N  and  not  z  or  0  .  This 
expression  was  arrived  at  by  taking  0  to  be  the  true  value  0(> 
of  the  process  generating  the  actual  observations  z.  However, 
this  is  not  significant  since  a  formula  for  evaluating  the 
probability  density  p(z,0Q)  as  a  function  of  both  z  and  0 
must  yield  the  correct  result  for  any  possible  z  whether  or 
not  it  in  fact  came  from  a  model  with  0q  the  true  parameter 
value . 


To  maximize  the  likelihood  function,  gradients  of  the 
log  likelihood  function  are  needed.  This  is  obtained  by  dif¬ 
ferentiating  the  log  likelihood  function  with  respect  to  each 
component  0j  of  0 
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3  In  p(z  ,6  ) 

36/ 

J 


1 

2 


N  i 

£  ltr(B.1 
i  =  l  1 


3  B  . 


36  . 
J 


i> 


rp  ,  3  B  ■ 

vT  B"1  ~  BT1  V  • 
-1  l  36 .  1  -l 


T'  1  3v  • 
+  2  b;  ^ 


(4.1-6) 


The  sensitivities  3B../30j  and  3v^/36^.  are  expressed  for  a  given 

6  .  in  terms  of  the  sensitivities  3F.  ,./30  .  and  3x.  ,./30  .  of 
J  t/i  J  -i/i  j 

the  state  estimation  error  covariance  and  state  estimate  by 
differentiating  Eq .  4.1-3  and  4.1-4.  This  gives 


3v  . 
-  i 

36  . 
J 


3H. 

l 


36j  -i/i - 1 


H. 

l 


3x  .  . .  , 
-l/i - 1 

36  . 

J 


(4.1-7) 


3  B  . 

_ l 

36  . 
J 


3H 


3  P 


P.,.  ,  HT  +  H.  — HlT  +  H . 

36  .  l/i-l  l  l  36 j  I  ] 


3H 


-Z  + 


3H  . 


i/i -1  36  .  '  36 

J  .1 


(4.1-8) 


The  sensitivities  of  the  state  estimate  and  its  covariance 

matrix  needed  in  the  above  equations  are  obtained  in  Section 

4.1.3  by  straightforward  but  lengthy  differentiation  of  all  the 

Kalman  filter  equations  with  respect  to  6  .  This  results  in  .j 

set  of  equations  similar  in  form  to  the  Kalman  filter  (and  in 

order  of  computation)  except  that  instead  of  propagating  the 

state  estimate  x.  and  its  covariance  matrix  P.  ,  the  sensi- 
-1/1  i/i 

tivities  3x^-/90j  and  3P^^/36j  are  propagated. 

Finally,  an  approximation  to  the  Fisher  information 
matrix  is  obtained  from  the  gradient  calculations.  The  (i.j) 
element  of  the  Fisher  information  is  approximated  as  (Kef.  20) 
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Propagation  of  sensitivities  of  the  state  estimate 
3x .  ,./30  .  and  its  error  covariance  matrix  3P.  ,./90  •  with  re- 
spect  to  a  parameter  0^  are  detailed  in  this  section.  This 
involves  straightforward  differentiation  of  the  Kalman  filter 
equations  detailed  in  Section  2.  First  the  state  sensitivites 
are  discussed  and  then  the  covariance  sensitivities  given. 


To  begin  with,  we  note  that  the  nonlinear  differential 
equations  (2.A.1)  are  parameter  dependent  so  that  we  write 
Fc(xc,0,t)  as  dependent  upon  the  parameters  0.  To  avoid  pos¬ 
sible  confusion  in  differentiation,  all  derivatives  of  F  with 

-c 

respect  to  xc  are  considered  with  xc  as  independent  of  0  the 

chain  rule  is  used  to  take  into  account  subsequent  dependence 

of  x  on  0 . 

-c  - 


The  discrete  and  continuous  state  estimates  x.  and  x 

-d  -c 

are  propagated  according  to  (2.5-3)  and  (2.A-5)  with  the  term 
deleted.  Straightforward  differentiation  with  respect  to 
0 j ,  noting  that  A  and  B  are  not  functions  of  0  .  ,  gives 


3 X  ( t  +A )  3x  ( t  )  3 h 

-Cx  o  _  -cx  o'  3D  h  +  n  — - 

30  .  *  30  .  30  .  -  30  . 

J  J  J  J 


( A  .  1  -  1  0  ) 
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3  x  .  [  ( vk )  °  ]  3x  ,[  ( v  k )  ]  3x  |(vk)  ] 

u _  -  a  _ - _  +  R  _ - _ 

ae  .  ae .  ae . 

j  j  j 


(4.1-11  ) 


where 


h  =  Ic^Ec^o*®)  ’  ®  •  t0) 


(4.1-12) 


a£c 

F  =  - — 

c  Q  T 

aic 


(ac(tG).  tQ) 


(4.1-13) 


and  D  is  the  matrix  function  (exp(FcA)  -  l)Fc  . 
atives  3h/80j  are  computed  by  the  chain  rule 


Partial  deriv- 


3F  3x  3  F 

3-/3eJ  =  SIT  +  5eT 

c  (E0.  0.  t) 


(4.1-14) 


Denoting  F  =  Fc  ,  the  computation  of  <t>(t,t^_^)  and 
D(t,t^  | )  can  be  done  along  with  their  partial  derivatives  as 


<b  (At ) 


“>  /  r' .  _  \  n 

*  ♦  E,  ^ 


(4.1-14) 


<*  /  r-A  .  \n 

»<“>  =  •  *  ?.  r^rb 


( 4 . 1  -  1  <i ) 


a<t>(At ) 

ae  . 

j 


y*  1  3 ( FA  t  ) n 

4-  n  !  ae  . 
n=l  J 


(4.1-17) 


with  a(FAt)  /ae.  computed  recursively  as 

3  (  FAt  )n  _  _  3  F  ,  r' »  .  a  n-  1  .  a(FAt)n_1 

—ae: - At  ae~.  (FAt)  4  FAt - acT - 

j  j  j 


(4.1-14) 
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To  accelerate  convergence  of  these  series,  4>  and  D  can  be  com¬ 
puted  using  these  series  for  a  time  increment  At  2~n  and  then 
4>(At)  and  D(At)  computed  recursively  by 


4>  (A  t  2"k)  =  [4>  (At  2"k_1 )  l2 


(4.1-20) 


D(At  2'k)  =  F_  1  ( I  -  4>  ( A  t  2"k_1))  (I  +  4>  ( A  t  2'k"1)) 


=  D(At  2_k_1)  (1  +  4>  ( A  t  2'k_1)) 


(4.1-21  ) 


.  1/  _  1/ 

and  the  partial  derivatives  34>(At  2  )/30j  and  3D(At  2  )/30  . 

computed  by  differentiating  these  two  recursive  relations. 


Sensitivities  of  the  state  estimation  error  covari¬ 
ance  matrix  can  be  conveniently  propagated  along  with  the 
state  sensitivities.  The  partials  of  the  state  covariance 
matrices  3Pc/30 ^  ^  ,  3Pccj/30j  are  obtained  by  differ¬ 

entiating  items  9  and  12  of  the  Filtering  Algorithm  in  Section 
2.5  as 


r  *  (Hr)  rA  *  ♦c  (§r)  *1 4  *Aw)  * 

J  J  J  J 


(4.1-22) 


3P  ,34>  3P 

_£d  _  p  *  _cd 

36.  36.  cd  c  36. 

J  J  J 


(4.  1-2.1) 


for  time  transition  (v(k+l))  to  (v(k+l))  and 


£  • 

J  '  J'  '  j'  '  J 


*  w.  <BPc/’ 


*  sr.  ,BPr/|T  ♦  sr. 


(4. 1-24 ; 
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*T  *  pcd(Hr)T*(Sf) *  pc(H:)1 


(4.1-25) 


for  time  transition  (v(k+l))  to  (v(k+l))°. 


The  remaining  relations  concern  the  Kalman  update. 

The  Kalman  gain  sensitivity  is 

3K 

30. 

J 

(HPHT  +  R)'1  (4.1-26) 


3P 

30. 


H1  +  P 


3HT 

30. 

l 


(HPH  +  R) 


-1  3H_ 
30. 

l 


PH  +  H 


3P 

30. 

l 


H  +  HP 


(8:) 


3R 

30. 

i 


Sensitivities  of  the  updated  state  estimate  and  its  error  co- 
variance  matrix  are  given  by  differentiating  items  17  and  18 
of  the  Filtering  Algorithm  in  Section  2.5: 


3x 

3 x 

.  3K 

3  v 

i  +  K  86. 

30  . 

ar: 

+  307 

J 

J 

J 

J 

3  P 

3  P 

3  K 

HP  -  K  Hr 

30  . 

30  . 

30  . 

KH 


3  P 
30 


(4.1-27) 


(4.1-28) 


4.2  OPTIMIZATION 

One  of  the  greatest  challenges  in  parameter  identi¬ 
fication  has  been  the  solution  of  numerically  i 1 } -condi t ioned 
or  even  ill-posed  maximization  problems.  This  is  a  problem 
receiving  considerable  attention  in  1  east  -  squares  methods  and 
much  less  attention  in  the  maximum  likelihood  case.  The  de¬ 
sign  of  general  purpose  algorithms  which  work  reliably  and 
efficiently  is  necessary  for  the  general  application  of  maxi¬ 
mum  likelihood  theory.  This  is  particularly  relevant  for  the 
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identification  of  catchment  parameters  where  the  differential 
equations  are  time  varying  and  involve  saturation  effects  with 
structural  changes  in  the  equations  at  overflow  of  various 
reservoirs.  Because  of  these  complications,  a  general  and 
reliable  procedure  is  required. 

Described  below  is  an  approach  based  upon  a  new  small 
sample  theory  of  ident i f i abi 1 i ty  which  is  usable  even  when 
some  functions  of  the  parameters  are  non i den t i f i abl e ,  i.e., 
when  perturbations  of  some  combinations  of  the  parameters  0 
produce  no  change  in  the  likelihood  function.  In  this  approach, 
the  nonident i f iable  parameters  are  first  determined  by  inspect¬ 
ing  the  Fisher  information  matrix.  Then  a  Levenburg-Marquardt 
optimization  procedure  is  used  in  the  subspace  of  parameter 
space  that  is  orthogonal  to  the  nonident i f i abl e  parameters  to 
compute  an  identifiable  set  of  parameters. 

U . 2 . 1  Ident i f iabi 1 i t y  Theory 

The  introduction  of  parametric  statistical  inference 
concepts  and  methods  by  R.A.  Fisher  (Ref.  21)  was  one  of  his 
major  contributions  to  statistics.  Unfortunately,  a  major 
difficulty  in  ident i f i abi 1 i t y  problems  has  been  an  overemphasis 
by  many  researchers  on  inference  about  a  parameter  vector  0 
rather  than  inference  about  the  class  F  =  Jp(z.O)}  of  probabil¬ 
ity  densities  i ndexed  by  the  parameter  0.  Most  definitions  of 
i den t i f i ab i 1 i ty  concern  properties  of  the  resulting  parameter 
estimates  rather  than  intrinsic  properties  of  the  class  F  in¬ 
dexed  by  the  parameter  0.  Indeed,  some  definitions  require  a 
hypothetical  infinite  sample  and  define  ident i f i abi 1 i t y  in  terms 
of  asymptotic  convergence  of  the  estimates  to  the  true  parameter 
values. 


This  overemphasis  on  the  parameter  values  rather  than 
the  class  of  probability  densities  indexed  by  the  parameters 
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has  developed  despite  the  early  and  fundamental  contribution 
of  Ref.  22.  This  paper  explicitly  defines  i den t i f i abi 1 i ty  of 
the  class  F  if  unique  parameter  values  produce  unique  probabili¬ 
ty  densities.  This  formulation  of  the  problem  of  specifying 
probability  models  for  statistical  inference  includes  the  iden- 
tifiability  problem  -  i.e.,  whether  or  not  the  specified  models 
are  "observationally"  unique.  Later  developments  in  the  litera¬ 
ture  seem  to  have  largely  overlooked  this  basic  concept  except 
for  a  few  econometric  papers  (see  Ref.  23  and  cited  references). 

In  the  approach  selected  for  this  study,  the  properties 
of  the  parameterized  class  {p(z,0),  6  t  0}  are  the  central 
issue.  The  above  definition  of  iden t i f i abi 1 i ty  (Ref.  22)  as 
formulated  for  the  parametric  case  by  Ref.  23  is  adopted: 

Two  parameter  points  0^  and  0^  are  said  to  be  observa¬ 
tionally  equivalent  if  p(z,0j)  =  p(z,09)  with  probabil¬ 
ity  1.  A  parameter  point  0^  is  said  to  be  global  1 y 
identifiable  if  there  is  no  other  0  c  0  which  is  ob¬ 
servationally  equivalent.  A  parameter  point  0j  is 
said  to  be  locally  identifiable  if  there  exists  an 
open  neighborhood  of  0j  containing  no  other  0  in  0 
which  is  observat ional ly  equivalent. 


This  approach  exploits  the  equivalence  between  local 
iden t i f i abi 1 i ty  and  full  rank  of  the  Fisher  information  matrix 
as  in  Ref.  23.  To  extend  this  connection  much  more  generally, 
a  powerful  new  result  on  existence  of  identifiable  reparamet er- 
izations  is  used  (Ref.  24 )  . 

Reparamet  erizat ion  Theorem  -  If  the  Fisher  informal  ion 
matrix  Ffl  of  a  paramet erizat ion  0  of  the  likelihood  (unction 
has  consTant  rank  h  in  a  neighborhood  of  a  point  0  of  parame¬ 
ter  space,  then  there  exists  a  reparamet  erizat  ion  tji(O)  such 
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that  (i|>  j  ,  .  .  .  )  is  locally  identifiable  and  the  likelihood 

function  is  not  a  function  of  +  ^  ,  .  .  .  .  Furthermore  the 

gradient  vectors 

(34>j/ae1 . 94ij/aek)T  for  j  =  h+l,...,k 

span  the  null  space  of  F„  . 

By  using  reparameterizat ions  resulting  in  a  nonsingu¬ 
lar  Fisher  information  matrix,  a  complete  characterizat i on  of 
local  ident i f iabi 1 i ty  by  the  Fisher  information  is  possible  in 
the  singular  case.  Previous  results  have  carefully  avoided 
reparameterizat ions  where  existence  is  not  trivially  guaran¬ 
teed  by  constraints,  etc.  Such  a  reparameterization  seems  to 
be  necessary  to  obtain  these  general  results. 

Another  aspect  of  the  ident i f iabi 1 i ty  approach  is  to 
exploit  the  special  structure  involving  the  Fisher  information 
matrix  to  devise  efficient  and  numerically  well  conditioned 
methods  for  maximizing  the  likelihood  function.  Using  the 
general  results  on  reparameterizations ,  it  is  possible  to  gen¬ 
eralize  and  to  make  precise  a  procedure  for  using  generalized 
inverses  in  the  method  of  scoring  when  the  Fisher  information 
matrix  is  singular  (Ref.  25).  Specifically  how  the  reparame- 
terization  result  is  useful  in  studying  the  special  structure 
of  the  maximum  likelihood  optimization  problem  is  discussed  in 
detail  in  the  following  section. 

4.2.2  Maximization  of  likelihood  Functions 

Lack  of  uniqueness,  i.e,  nonident i f iabi 1 i ty ,  manifests 
itself  as  i 1 l -condi t ioni ng  in  computation  of  least  squares  or 
maximum  likelihood  estimates.  Even  when  the  parameters  un¬ 
identifiable,  i 1 1 -condi t ioni ng  often  arises  because  of  the 
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considerable  difference  in  sensitivity  of  the  likelihood  or 
squared  error  functions  to  changes  in  different  parameters. 

Some  past  work  used  the  special  structure  of  least  squares 
problems  to  devise  efficient  optimization  methods  which  recog¬ 
nize  any  nonuniqueness  or  i 1 lcondi t ioning  and  solve  for  the 
parameters  within  the  equivalence  due  to  the  nonuniqueness 
(Refs.  26,  27).  Very  little  progress  along  these  lines  has 
been  made  in  the  maximum  likelihood  problem  although  a  very 
general  method  has  been  proposed  (Ref.  25).  The  maximum  like¬ 
lihood  method  provides  practical  parameter  estimation  and  test 
of  hypothesis  procedures  for  many  complex  random  processes, 
and  also  provides  the  needed  approximate  distribution  theory 
which  is  not  generally  available  with  alternative  procedures. 

Most  maximum  likelihood  methods  require  that  the  class  of  mod¬ 
els  be  reparameterized  uniquely  so  all  parameters  are  identifi¬ 
able.  There  are  no  conditions  given  for  when  such  a  reparame¬ 
terization  is  possible.  To  answer  questions  of  existence  and 
to  actually  reparameterize  involves  solving  a  system  of  nonlinear 
partial  differential  equations. 

The  method  proposed  in  Ref.  25  presumes  that  there 
exists  a  reparameterization  for  which  the  Fisher  information 
matrix  is  nonsingular  and  evidently  equates  such  nonsingularity 
to  ident i f iabi 1 i ty  although  any  rigorous  discussion  or  even 
definition  of  ident i f iabi 1 i ty  in  the  nonlinear  case  is  lacking. 

It  is  argued  that  the  method  of  scoring  (using  the  Fisher  infor¬ 
mation  matrix  in  place  of  the  Hessian  in  a  Newton  type  algorithm) 
can  be  implemented  entirely  in  the  original  nonident i f iable 
parameterization  using  the  pseudoinverse  of  the  Fisher  infor¬ 
mation  matrix  to  restrict  the  maximization  to  a  locally  identi¬ 
fiable  subspace  of  parameter  space.  This  would  avoid  any  need 
to  reparameterize  in  terms  of  an  identifiable  set  of  parameters. 
Presently  the  only  alternative  method  to  preclude  i dent i f i abi 1 i t y 
difficulties  is  to  actually  carry  out  such  a  reparamet  eriz.it  ion 
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which  is  often  exceedingly  difficult  if  at  all  feasible  computa¬ 
tionally.  Because  of  the  complexity  of  the  catchment  model 
discussed  above,  determining  such  a  reparameterization  almost 
certainly  would  not  be  feasible. 

One  of  the  objectives  of  the  present  approach  is  to 
use  the  special  structure  of  the  maximum  likelihood  estimation 
problem  to  devise  efficient,  numerically  accurate  and  stable 
maximization  methods.  In  particular  the  method  of  Kef.  2r>  can 
be  shown  to  work  very  generally  utilizing  the  new  result  de¬ 
scribed  in  Section  4.2.1  which  guarantees  the  local  existence 
of  an  identifiable  reparameterizat ion  whenever  the  Fisher  in¬ 
formation  matrix  has  constant  rank  locally.  This  also  general¬ 
izes  the  equivalence  of  ident i fiabi 1 i ty  and  full  rank  of  the 
Fisher  information  matrix  (Ref.  23)  to  the  reduced  r  -k  case 
(the  null  space  is  the  local  equivalence  class). 

The  method  of  scoring  is  not  only  attractive  in  re¬ 
moving  ident i fiabi 1 i ty  problems,  but  it  has  several  other  at¬ 
tractive  computational  features.  In  a  number  of  problems  the 
Fisher  information  can  be  computed  from  the  gradient  computa¬ 
tions  with  little  additional  work.  Thus  Hessian  information 
is  obtained  from  the  gradients  without  additional  computation. 

It  has  been  found  in  practice  that  this  Hessian  approximation 
gives  excellent  approximation  to  the  eigenvectors  which  domi¬ 
nate  the  numerical  behavior  of  approximate  Hessian  methods. 

This  suggests  further  special  structure  of  the  problem  since 
Hessian  approximations  in  i 1 lcondi t ioned  cases  usually  result 
in  poor  algorithm  performance  unless  there  is  some  special 
st  ructure . 

4.2.3  Quadratic  Algorithm  for  i  liable  Parameters 

The  above  theory  provides  a  basis  for  an  efficient 
and  well  behaved  algorithm  even  though  there  may  be  non i dent i- 
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fiable  parameters  present.  The  first  step  is  to  perform  an 
eigenvalue-eigenvector  decomposition  of  the  Fisher  information 
matrix  F  in  the  form 


F  =  U  A  UT  (4.2-1  ) 

T 

where  the  columns  of  U  are  orthogonal  vectors  so  U  U  =  I  and 

where  A  is  diagonal  with  elements  arranged  in  descending  order 

Aj  >  ...  ^A^^O.  The  algorithm  should  automat i cal ly  check 

that  F  is  positive  definite  to  detect  numerical  difficulties 

that  might  have  occured  in  computing  F.  To  determine  identi- 

- 1 2 

fiability,  a  threshold  c,  say  10  for  double  precision  with 
15  decimal  places,  is  set  and  any  eigenvalues  A^  £  eAj  are 
mod i f i ed  as 


(  A .  if  A  .  >  e A . 

A.  .  =  1  1  1  (4.2-2) 

1  mod  0  if  A.  <  eA. 

'  l-l 

The  eigenvectors  corresponding  to  the  modified  eigenvalues  set 
to  zero  specify  the  locally  nonident i f i abl e  linear  combinations 
u .0  of  the  parameters.  The  optimization  is  constrained  to  the 
subspace  of  linear  combinations  of  parameters  that  are  ortho¬ 
gonal  to  these  non i den t i f i abl e  combinations,  i.e.,  constrained 
to  the  subspace  of  parameter  space  spanned  by  the  eigenvectors 
Uj...,u^.  corresponding  to  the  unmodified  eigenvalues. 

A  quadratic  algorithm  is  one  which  uses  local  quadrat¬ 
ic  information  about  the  log  likelihood  function,  i.e.,  the 
gradient  (first  partial  der ivat ives )  and  the  Hessian  (second 
partial  derivatives).  The  identi f iabi 1 i ty  theory  guarantees 
that  the  null  space  of  the  Hessian  will  exactly  coincide  with 
that  of  the  Fisher  information  matrix  so  that  perturbations  ot 
linear  combinations  in  the  null  space  of  the  Fisher  information 
matrix  wi 1 1  have  no  effect  upon  the  likelihood  (unction.  Since 
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the  log  likelihood  function  is  in  general  only  locally  quadratic 
some  procedure  is  needed  to  constrain  the  step  A0  in  parameter 
space  to  a  region  in  which  this  quadratic  approximation  is 
a  reasonable  one.  An  excellent  procedure  is  the  so-called 
Levenburg-Marquardt  procedure  (Ref.  28)  which  is  discussed  in 
a  statistical  setting  in  Ref.  29  as  an  optimal  step  that  maxi¬ 
mizes  the  quadratic  approximation  subject  to  a  fixed  step  length 
If,  in  addition,  the  optimization  is  constrained  to  be  orthog¬ 
onal  to  the  nonident i f iable  parameters,  then  the  algorithm  has 
the  form 


“  =  -(Fmod  +  Mmod>  l  ln  P 


(4.2-3) 


where 


V  In  p  =  vector  of  gradients  8  In  p(z,0)/30^ 

F  ,  =  U  A  ,  UT 
mod  mod 


^od  =  U  fi  U 


tT 


1  if  \mod  >  0 


fi  is  diagonal  with  w .  .  =  . 

11  0  if  A .  ,  =  0 

'  imod 

A  is  a  step  length  parameter 
t  is  the  pseudoinverse  operation 

The  calculation  is  implemented  using  only  vector  products  as 


k  L(t  In  p) 
A0  =  -  £  1 

i  =  1  i 


A.+A  -i 


(A.2-A) 


where  V  In  p  is  the  gradient  vector  9  In  p/90. 


Trial  steps  are  made  with  an  initial  value  of  A  from 
initialization  of  tKe  program  or  from  the  last  iteration.  For 


A  -  1  8 
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the  trial  step  0^  +  j  =  6^  +  40 ,  the  log  likelihood  function  is 
evaluated  and  the  quadratic  prediction  criterion 


=  1  - 


>redicted  lo^ 
actual  log 


likelihood 

likelihood 


(A.  2-5) 


is  calculated  where  the  predicted  log  likelihood  In  p  is  based 
upon  the  quadratic  approximation 

In  p(e..,)  =  In  p(0.)  +  48T  V  In  p  +  48T  F  .40 
—l  +  l  —  i  —  —  -  mod  — 

(4.2-6) 

The  step  length  parameter  A  is  adjusted  by  powers  of  2  until 

.25  <  p  <  .75  (4.2-7) 

or  until  A  <  A^/8  in  which  case  A  is  set  equal  to  zero  and  a 
modified  Newton  step  results. 

The  above  optimization  procedure  is  implemented  in 
TM 

the  PARAIDE  computer  software.  This  procedure  has  been  used 
on  a  number  of  very  complex  nonlinear  parameter  identification 
problems  in  the  past  with  very  good  numerical  behavior  and 
convergence . 


The  rate  of  convergence  of  the  quadratic  algorithm 
depends  upon  the  goodness  of  the  quadratic  approximation  to 
the  log  likelihood  function.  This  quadratic  approximation  can 
be  improved  by  choice  of  an  appropriate  parameterization  of 
the  likelihood  function  which  stabilizes  the  Fisher  information 
matrix,  that  is,  which  results  in  a  nearly  constant  matrix 
with  changes  in  the  parameters  (Ref.  30).  Such  a  stabilizing 
parameterization  also  improves  the  approximate  distribution 
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theory  on  the  parameter  estimates  and  results  in  very  nearly 
normal  distributions  of  the  estimates  even  for  quite  small 
samples  (Ref.  31).  This  is  true  in  spite  of  the  very  nonline¬ 
ar  nature  of  the  catchment  model  so  that  with  an  appropriate 
reparameterization  the  problem  is  nearly  quadratic. 

Stabilizing  reparameterizations  are  obtained  from  a 

given  parameterization  6  and  its  Fisher  information  matrix  FQ 

—  u 

by  finding  a  reparameterizat ion  0(6)  of  6  with  a  gradient  A^q 
such  that 


q,  =  <v>'1T  w’’1 

is  nearly  constant  independent  of  the  point  0  of  parameter 

space,  where  (J>  is  the  matrix  of  partial  derivatives  with 

(i,j)  element  8i {k/30 j  .  This  is  in  general  difficult  to  do, 

but  in  many  cases  there  are  simple  reparamet erizat ions  which 

will  yield  considerable  improvement.  Two  examples  are  the 

2 

variance  parameter  o  and  the  correlation  coefficient  p  which 
are  improved  by  the  reparameterizations 

0(a2)  =  In  o2  (6.2-9) 

<(>(p)  =  ^  In  =  arc  tanh  (p)  (6.2-10; 


In  general  for  a  single  parameter  6  the  stabilizing  reparame¬ 
terization  is  given  by  integrating  Eq .  6.2-8  as 


0(6) 


(Fq)1/2  d6 


6.2.5  Statistical  Convergence  Criterion 


(6.2-11  ) 


A  convergence  criterion  is  used  to  decide  when  the 
maximum  of  the  likelihood  function  has  been  found  with  suf- 
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ficient  precision.  The  objective  is  to  make  the  error  0k-0 
between  the  computed  value  0^  at  iteration  k  and  the  true  maxi¬ 
mum  0  small  relative  to  the  expected  error  in  0  due  to  sampling 
variability.  Thus  the  convergence  measure 


(eR-0)T  [  cov  ( 0  )  ]  ~ 1  <ek-0)  <  (6k-6)  F  (0  )  (0  k~0  )  =  6 

(A. 2-12) 

is  used,  where  the  inequality  follows  from  the  Cramer-Rao  lower 
bound  (Ref.  32).  Since  the  function  is  usually  nearly  quadratic 
close  to  the  solution,  we  have 

i  *  2k  -  =  -  Fk  tk(ln  p>  ''■•2-13) 


and  using  the  notation  VR  =  V  R(ln  P)  yields 


6  $  il  Fl  ?k 


(A.2-1A  ) 


Also  the  error 
In  p(0)  is 

£n  p ( 0 )  - 


in  computing  the  value  of  the  function 


£n  p(ek)  =  vk(0-0k)  +  (e-ek)T  H(e-ek) 

5  +  <?k+re-k>T  H<6-k+i-?k> 

=  tk<«k,i-&k>  +  <WVT  "'“k.i-'k1 


T  t  T  t  t 

=  YkF TVk  +  V*FTHFTVk 


(  A  .  2  -  1 1 ) 


where  H  denotes  the  Hessian  matrix.  Thus  on  the  average,  the 
maximum  value  £n  p(0)  is  underestimated  by 


E | £n  p(0)  -  £  n  p<  0 k )  )  =  2  VR  Ft  yk  =  26 


(  A  .  2  -  1  (> ) 


A  -  2 1 
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The  quantity  -2£n  p(0)  is  basic  in  tests  of  hypotheses.  It 
can  be  shown  that  the  addition  of  each  unnecessary  parameter 

to  the  model  will  decrease  the  value  of  the  test  statistic 

'  2 

-2£n  p(0)  by  a  chi-square  random  variable  x-i  on  one  degree  of 

2  ^ 

freedom  (Ref.  16).  Since  Var(xt)  =  2,  the  error  in  computing 

i  p 

-2  £n  p ( 6  )  as  -2  £n  p(0^)  is  46 .  Thus  for  6  =  10  ,  the  error 
in  computing  2  £n  p(6)  is  small  compared  with  the  expected 
sampling  variation. 


4.3  RESULTS 

The  maximum  likelihood  parameter  estimation  methods 
are  demonstrated  on  precipitation  and  discharge  measurements 
simulated  using  a  NWS  supplied  catchment  model.  The  parameter 
identifiability  question  is  discussed,  and  then  the  local  be¬ 
havior  of  the  algorithms  and  nature  of  the  parameter  estimates 
described . 

4.3.1  Parameter  Identifiability 

The  identifiability  of  catchment  model  parameters  is 
illustrated  in  this  section  for  two  contrasting  cases  of  heavy 
and  moderate  rainfall.  The  rainfall  and  channel  discharge 
data  was  simulated  for  Bird  Creek  basin  over  a  one-month  period. 
The  hydrologic  structures  excited  by  these  rainfall  records 
are  compared  in  Table  4.3-1.  The  rainfall  records  are  illus¬ 
trated  in  Figs.  2.2-4  and  2.6-3  for  the  moderate  and  heavy 
rainfall  cases  respectively.  The  various  catchment  model 
state  histories  for  the  heavy  rainfall  case  are  illustrated  in 
Figs.  2.6-4  through  2.6-19. 
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TABLE  A. 3-1 

COMPARISON  OF  MODERATE  AND  HEAVY  RAINFALL  CASES 


CASE 

CATCHMENT^\_ 

ELEMENT 

MODERATE 

RAINFALL 

HEAVY 

RAINFALL 

Maximum  upper-zone 
tension  water  content 

3  time 

6  times 

Surface  runoff 

none 

3  times 

Maximum  lower- zone 
tension  water  content 

1  time 

none 

Upper-zone  free-water 
content 

very 

little 

nonzero 
half  of  time 

The  ident i f iabi 1 i ty  of  parameters  for  the  heavy  rain¬ 
fall  case  is  shown  in  Table  4.3-2.  The  first  two  columns  show 
the  catchment  model  parameters  and  their  true  values  used  in 
the  simulation  of  channel  discharge  data.  The  approximate 
standard  deviations  and  correlations  of  the  parameter  estima¬ 
tion  errors  were  calculated  from  the  inverse  Fisher  information 
matrix.  Wh>jn  all  fifteen  parameters  are  simultaneously  esti¬ 
mated,  only  x”  and  x^  have  standard  deviations  around  1  percent 
while  the  parameters  x.^,  d^,  d^  ,  d£  ,  p  ,  and  a^  are  less  than 
about  10  percent,  and  the  remaining  parameters  have  errors  ol 
the  order  of  the  parameter  values  themselves.  If  however  all 
parameters  but  one  were  known  or  presumed  known,  then  the  stan¬ 
dard  deviation  due  to  sampling  variability  in  estimating  only 
that  one  parameter  is  reduced  as  shown  in  column  4.  In  some 
cases  the  reduction  is  nearly  two  orders  of  magnitude  and  gen¬ 
erally  involves  high  correlations  with  other  parameters  when 
estimating  all  fifteen  parameters.  These  high  correlations 
indicate  near  deterministic  dependence  between  the  parameters 
inv'olv'ed  and  implies  that  corresponding  changes  in  these  param¬ 
eters  have  almost  no  effect  upon  the  predicted  observations. 
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TABLE  4. 3-2 

IDENTI F I ABILITY  OF  CATCHMENT  PARAMETERS 
FOR  HEAVY  RAINFALL  CAST 


PARAMETER 


TRUE 

VALUE 


STANDARD  DEVIATION 

15  PARAMETERS 
ESTIMATED 

1  PARAMETER 
ESTIMATED 

1.08 

0.984 

.  188 

0.019 

9.8 

1.20 

57 

0.88 

9. 1 

0.  19 

3.9  E-4 

5.1  E-5 

5.8  E-5 

5.0  E-6 

1.1  E-3 

7.0  E-5 

42 

1.4 

1.2 

0.022 

0.009 

0.0078 

0.33 

0.043 

0.015 

0.0046 

0.00094 

0.00049 

0.88 

0.30 

HIGH 

CORRELATIONS, 
15  PARAMETERS 
ESTIMATED 


0.87  (x“  x?) 

.5  4 


0.88  (x°,  v) 


-0.90  (x°,  d'p 


0.92  (x4,  a) 


0.993  (y,  a) 


120  mm 
15  mm 
160  mm 
140  mm 
14  mm 

1.486  E-2  1/hr 
5.452  E-4  1/hr 
5.612  E-3  1/hr 
48 
2.1 
0.02 
3.55 
0.17 
0.001 
0.0 


In  the  case  of  estimating  only  one  parameter  with  all  others 
known,  the  error  is  less  than  about  1  percent  in  all  parameters 
except  pp  aj,  and  H2~  The  more  difficult  to  determine  param¬ 
eters  are  associated  with  the  percolation  function  (2.2-18)  in 
both  cases  of  estimating  one  parameter  alone  or  a  1 1  fifteen 
s imu 1 1  aneous 1 y . 
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The  ident i f iabi 1 i ty  of  parameters  for  the  moderate 
rainfall  case  is  shown  in  Table  4.3-3.  It  is  qualitat ively 
similar  to  that  of  the  heavy  rainfall  case  except  that  the 
estimation  error  standard  deviations  are  larger,  by  orders  of 
magnitude  in  some  instances,  and  the  high  correlations  are 
more  numerous  and  some  very  high  (greater  than  0.99).  When 
all  fifteen  parameters  are  estimated,  only  x°  can  be  deter¬ 
mined  to  about  1  percent,  X£  to  about  25  percent,  and  the  rest 
are  not  identifiable.  Much  of  this  difficulty  is  due  to  the 
inability  to  differentiate  between  the  effects  of  perturbations 
in  different  parameters  upon  the  discharge  measurement.  This 
is  apparent  since  individual  estimation  of  a  parameter  with 
the  others  known  results  in  an  error  of  less  than  several  per¬ 
cent  except  for  the  parameters  p^  and  Note  as  before  that 

the  most  dramatic  reduction  in  estimation  error  between  the 
cases  of  estimating  only  one  parameter  and  estimating  all  fif¬ 
teen  simul taneously  occurs  when  very  high  correlations  are 
involved . 


The  conclusions  apparent  from  the  ident i f iabi 1 i t y 
results  for  the  moderate  and  heavy  rainfall  cases  above  are 
that : 

•  Sufficient  rainfall  to  excite  all  basin 
dynamics  is  required  if  parameters  are 
to  be  estimated  simultaneously 

•  Some  of  the  parameters  associated  with 
the  percolation  function  are  not  identi¬ 
fiable  even  if  all  other  parameters  were 
known  exactly. 


The  i den t i f i ab i 1 i t y  of  parameters  for  a  given  basin  will  thus 
depend  very  strongly  upon  the  available  records  ot  storms  or 
high  snowmelt.  Just  increasing  the  length  of  the  data  used 
will  not  significantly  improve  the  i den t i f i ab i 1 i t y .  What  is 
needed  is  more  data  involving  the  excitation  of  the  basic 
dynamics  in  different  ways. 
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TABLE  4.3-3 

I DENTL FI ABILITY  OF  CATCHMENT  PARAMETERS 
FOR  MODERATE  RAINFALL  CASE 


PARAMETER 


TRUE 

VALUE 


120  mm 
15  mm 
160  mm 
140  mm 
14  mm 

1.486  E-2  1/hr 
5.452  E-4  1/hr 
5.612  E-3  1/hr 
48 
2.  1 
0.02 
3.55 
0.17 
0.001 
0.0 


STANDARD  DEVIATION 

15  PARAMETERS 
ESTIMATED 

1  PARAMETER 
ESTIMATED 

'  "  ~  '  ■  ' 

1.4 

1  .0 

3.5 

0.56 

3600. 

2.6 

3200. 

2.  1 

120. 

0.72 

14.8  E-2 

0.077  E-2 

7.6  E-4 

0.06  E-4 

5.7  E-3 

0.21  E-3 

2200. 

3.1 

56. 

0.042 

0.42 

0.063 

50. 

0.039 

7.9 

0.007 

0.0025 

0.0009 

8.  7 

0.22 

HIGH 

CORRELATIONS, 
15  PARAMETERS 
ESTIMATED 


0.8b  (x°,  x£ ) 


0.84  (x  ,  (1  ) 
3  u 

-0.87  (x®,  x?) 

J  4 

-0.84  (x!J,  dg) 
0.83  (x“,  d^l 
-0.85  (x(3\  p) 
0.84  (x3  ,  a 1  ) 
-0.86  (u ,  a  j ) 
0.95  (y,  iv) 


-0.998  (d  ,  (.i) 
u 

0.88  (d’,  p) 

-0.998  (p,  a  1  I 

0.9999  (a,,  d  I 
I  u 

-0.98  (x®,  s) 


4.3.2  Demons  t ra t ion  of  Al go r i t  hm 


The  local  behavior  of  the  parameter  identification 
algorithm  near  the  maximum  of  the  likelihood  function  is  de¬ 
scribed.  The  global  behavior  from  various  initial  parameter 
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estimates  will  require  more  extensive  investigation  but  is 
expected  to  reproduce  the  generally  robust  and  efficient  be¬ 
havior  of  quadratic  hill  climbing  methods  (Ref.  33). 

The  case  considered  is  the  heavy  rainfall  case  dis¬ 
cussed  in  Section  4 . 3 . 1  and  the  precipitation  record  is  shown 
in  Fig.  2.6-3.  To  reduce  the  computer  time  required,  only  the 
first  10  days  were  used  and  the  three  parameters  x ° ,  x°,  and 
d^  were  simultaneously  estimated.  The  computer  CPU  time  on  an 
IBM  370/3031  using  a  PL/1  optimizing  compiler  was  about  3  min 
per  iteration  on  the  parameter  values  and  was  approximately 
the  number  of  parameters  estimated  (3  in  this  case)  times  the 
CPU  time  for  a  Kalman  filter  run.  There  is  great  potential 
for  considerably  reducing  the  required  CPU  time  by  exploiting 
the  special  problem  structure  especially  in  the  initial  stages 
where  the  parameter  estimates  are  far  away  from  the  maximum 
likelihood  estimates.  When  the  maximum  is  approached,  the  log 
likelihood  is  approximately  quadratic  and  convergence  acceler¬ 
ates  very  rapidly  as  discussed  below-.  Thus  few  iterations  are 
needed  near  the  solution. 

Table  4. 3-4  gives  the  result  of  four  iterations  of 

the  algorithm  starting  from  the  true  parameter  values.  These 

starting  values  are  about  17,.,  157...  and  107.  respectively  from 

the  maximum  likelihood  for  xV ,  x1' ,  and  d  .  The  standard  devi- 

12  u 

at  ions  show  little  change  for  Xj  but  moderate  change  for  x,; 
and  d^.  This  suggests  that  the  errors  in  estimating  x^  and  d^ 
have  a  probability  distribution  slightly  different  from  normal. 
A  reparameterizat ion  may  he1  ;  this  problem. 

The  changes  in  the  likelihood  function  value  are  due 
almost  entirely  to  changes  in  the  quadratic  term  which  is  the 
suin  of  normalized  innovations.  If  the-  parameter  est  i mates 
were  noramlly  distributed,  then  twice  the  difference  of  the 
log  likelihood  functions  evaluated  respectively  at  the  true 
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values  and  maximum  likelihood  estimated  values  of  the  parameters 
would  be  a  chi-square  random  variable  on  3  degrees  of  freedom 
(Ref.  32).  This  difference  is  12.45  which  is  statistically 
significant  at  the  0.01  level,  i.e.,  a  departure  this  large 
would  occur  less  than  1  percent  of  the  time  under  the  normal¬ 
ity  assumption.  This  indicates  significant  non-normality  and/ 
or  a  nonzero  mean  due  to  nonlinearities  in  the  catchment  and 

~  T 

rainfall  models.  Also,  the  quantity  (6-0)  F(0-0)  would  be  a 
chi-square  random  variable  on  3  degrees  of  freedom  if  the  param¬ 
eter  estimates  were  distributed  normally  with  zero  mean.  This 
quantity  is  18.17  and  9.68  respectively  for  the  Fisher  informa¬ 
tion  matrix  evaluated  at  the  true  parameter  0  and  at  the  maximum 
likelihood  estimates  0.  This  further  indicates  some  departure 
from  the  theoretical  distribution  of  the  estimates.  These 
quantities  would  be  expected  to  be  smaller  than  about  7.8  ninety- 
five  percent  of  the  time.  Such  a  departure  is  easily  caused 
by  a  standard  deviation  being  wrong  by  a  factor  of  2.  The 
above  discussion  indicates  that  the  computed  sampling  distribu¬ 
tion  as  indicated  by  the  Fisher  information  and  standard  devia¬ 
tions  are  accurate  within  a  factor  of  2. 

Further  nonlinear  effects  are  apparent  in  the  discharge 
and  precipitation  innovations.  The  precipitation  innovations 
in  Fig.  4.3-1  show  large  departures  around  4  days.  These  large 
departures  are  due  to  the  nonlinear,  nongaussian  precipitation 
model  whose  probability  density  has  heavy  tails  so  that  devia¬ 
tions  of  5  standard  deviations  are  not  too  unlikely.  This 
accounts  for  the  two  large  peaks  in  Fig.  4.3-4  which  displays 
the  innovations  normalized  by  their  inverse  covariance  matrix. 

Due  to  nonlinearities  in  the  precipitation  and  catchment  models, 
a  deterministic  filtering  error  appears  in  the  discharge  innova¬ 
tions  as  shown  in  Fig.  4.3-2  at  the  true  and  estimated  parameter 
values.  This  deterministic  function  is  somewhat  obscured  by 
the  innovations  noise  as  shown  in  Fig.  4.3-3,  however  in  other 
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TIME  (DAYS) 

Figure  4.3-1  Precipitation  Innovations 


discharge  innovations 


Time  (DAYS) 


Figure  4.3-2  Deterministic  Component  of  Channel  Discharge 
Innovations  for  True  (Solid  Line)  and  Ksti- 
mated  (Dashed  Line)  Parameter  Values 
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Figure  4.3 


DISCHARGE  INNOVATIONS 


3  Discharge  Innovations  for  True  (Solid  Line) 
and  Estimated  (Dashed  Line)  Parameter  Values 

NORMALIZED  INNOVATIONS 


Normalized  Innovations  for  True  (Solid  Line) 
ar.d  Estimated  (Dashed  Line)  Parameter  Values 


Figure  4.3-4 
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cases  it  can  be  expected  that,  relative  to  the  noise,  very 
large  deterministic  components  can  occur  in  the  discharge  in- 
novations  due  to  nonlinearities.  Except  for  the  peaks  due  to 
these  nonlinearities,  the  normalized  innovations  look  quite 
reasonable . 

The  convergence  behavior  of  the  algorithm  is  indi¬ 
cated  by  the  convergence  criterion  and  confirmed  by  the 
gradients  and  gradient  norm  (root  mean  square  of  gradients). 

The  convergence  criterion  exhibits  a  charac t eri s t i c  linear 

convergence  where  the  convergence  criterion  6  is  reduced  by 

_2 

approximately  the  same  factor  10  on  each  iteration.  Con¬ 
vergence  is  achieved  for  practical  purposes  on  the  second 

_2 

iteration  (6  <  10  ),  but  the  algorithm  was  continued  to 

demonstrate  the  character  of  the  convergence  and  precision 
of  the  computations. 

The  demonstration  of  the  maximum  likelihood  algorithm 
shows  the  potential  for  its  use  in  NWS  river  forecast  system. 
The  algorithm  behaved  efficiently  and  robustly  in  the  presence 
of  nonlinear  dynamics  and  nongaussian  noise.  Reasonable  values 
for  the  parameter  estimates  were  obtained  which  were  close  to 
the  error  predicted  by  the  approximate  sampling  distribution 
theory.  More  extensive  testing  of  the  algorithm  is  required 
to  determine  how  general  these  conclusions  are. 
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5. 


SUMMARY  AND  CONCLUSIONS 


Applications  of  modern  estimation  and  filtering  theory 
to  the  requirements  of  the  National  Weather  Service  (NWS)  were 
investigated  to  assess  their  potential  for  improving  river 
flow  forecasting  and  catchment  model  calibration. 

The  work  was  organized  into  three  principal  tasks: 

•  Issues  in  filter  design 

•  State-space  model  development  for  unit 
hydrographs 

•  Parameter  identification  for  catchment 
modeling. 


5.1  ISSUES  IN  FILTER  DESIGN 


The  applications  of  Kalman  filtering  techniques  to 
hydrologic  forecasting  required  the  development  of  catchment 
state-space  models.  These  were  obtained  by 


•  Developing  continuous-time  nonlinear 
state-space  equations  for  the  Sacramento 
Soil  Moisture  model 

•  Modifying  the  distribution  of  the  perco¬ 
lation  to  the  lower  zone  to  avoid  slight 
inconsistencies  in  the  Sacramento  model 

•  Modeling  channel  routing  with  reduced- 
order  state-space  models  for  unit  hydro¬ 
graphs 
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•  Using  state  augmentation  techniques  to 
combine  the  soil  moisture  and  channel 
routing  systems  to  generate  a  complete 
catchment  model. 

An  extended  Kalman  filter  for  the  estimation  of  the 
state  of  the  catchment  system  and  the  prediction  of  the  dis¬ 
charge  to  the  channel  was  designed  and  its  performance  tested 
with  simulated  data.  Excellent  agreement  between  true  and 
forecasted  flows  was  obtained  even  under  surface  runoff  condi¬ 
tions.  The  results  indicate  that  the  extended  Kalman  filter¬ 
ing  technique  constitutes  a  very  well  behaved  procedure  for 
the  practical  forecasting  of  the  discharge  from  a  basin. 


5.2  STATE-SPACE  MODEL  DEVELOPMENT  FOR  UNIT  HYDROGRAPHS 

The  canonical  variate  method  of  deriving  reduced- 
order  state-space  models  of  unit  hydrographs  gives 


•  Optimal  reduced-order  models  in  terms  of 
weighted  squared  prediction  error 

•  An  automatic  procedure  suitable  for  com¬ 
puter  implementation 

•  A  computationally  efficient  and  numerical¬ 
ly  stable  algorithm. 


The  modeling  of  NOAA/NWS  supplied  unit  hydrographs  indicates 
that 


•  The  sum  square  error  criterion  is  superior 
to  the  correlation  criterion  in  producing 
good  low  order  approximations  to  the 

unit  hydrograph 

•  Considerable  reduction  in  state  order 
was  typical  -  from  order  ten  or  twenty 
to  order  three  or  five 
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•  Multirate  unit  hydrographs  were  no  more 
difficult  to  approximate  than  those  with 
the  same  input-output  rate. 

This  reduced  order  state  space  modeling  effort  was  highly  suc¬ 
cessful  in  achieving  the  objectives  of  the  study.  One  inherent 
difficulty  in  multirate  unit  hydrographs  was  discovered  -  cer¬ 
tain  particular  unit  hydrographs  introduce  spurious  high  fre¬ 
quencies  in  the  output  as  a  result  of  a  low  frequency  sampled 
sine  wave  input.  The  need  for  multirate  unit  hydrographs  is 
completely  avoided  by  using  the  continuous  catchment  model  de¬ 
rived  in  this  report  and  integrating  to  obtain  a  discrete  time 
model  operating  at  the  same  rate  as  the  channel  routing  model. 


5.3  PARAMETER  IDENTIFICATION  FOR  CATCHMENT  MODELING 

The  application  of  maximum  likelihood  methods  to  the 
catchment  model  parameter  identification  problem  provides  an 
initial  evaluation  of  its  potential  use  by  NOAA/NWS.  The  maxi 
mum  likelihood  methods  described  in  this  study  provide 

•  A  determination  of  parameter  identifi- 
ability 

•  A  robust  optimization  algorithm  which  is 
immune  to  parameter  nonidenti f iabi 1 i ty 

•  Estimates  of  the  identifiable  parameters 
and  their  estimation  error  covariances 

•  An  automatic  procedure  suitable  for  com¬ 
puter  implementation. 

These  methods  were  used  on  simulated  precipitation  and  channel 
discharge  data  generated  using  a  NOAA/NWS  supplied  catchment 
model  for  the  Bird  Creek  basin.  Identification  of  the  catch¬ 
ment  model  parameters  indicates  that 
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•  Parameter  i dent i f iabi 1 i ty  is  determined 
by  the  extent  to  which  the  data  excites 
the  relevant  hydrologic  structures 

•  The  optimization  algorithm  converges 
very  rapidly  near  the  solution. 

The  maximum  likelihood  method  has  been  demonstrated  as  a  power¬ 
ful,  automatic  procedure  with  potential  for  wide  spread  use  in 
the  fitting  of  NOAA/NWS  catchment  models. 

5 . 4  RECOMMENDATIONS 


There  are  several  important  areas  for  future  investi¬ 
gation  suggested  by  this  study  in  the  identification  of  catch¬ 
ment  model  parameters 


•  Initialization  and  improved  computational 
efficiency  in  recursive  refinement  of 
parameter  estimates 

•  Tests  of  hypotheses  between  alternative 
hydrological  model  structures 

•  Use  of  robust  methods  on  data  for  han¬ 
dling  outliers 

•  Application  of  both  the  Kalman  filtering 
and  parameter  identification  to  a  number 
of  data  sets  for  a  variety  of  basins. 

Such  future  study  would  demonstrate  the  generality  with  which 
these  methods  apply  to  the  NOAA/NWS  operational  hydrologic 
forecasti ng . 
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APPENDIX  A 


The  constraint  on  the  ratio  of  free  to  tension  water 
for  the  upper  zone  is  analyzed  in  this  appendix.  It  is  shown 
below  that  this  constraint  is  superfluous  for  all  basins  examined. 

When  the  constraint  is  removed,  the  state-equat ions 
for  the  upper  zone  can  be  derived  as  follows  (Ref.  2).  Con¬ 
sider  first  the  tension-water  element.  Two  quantities  can 
affect  the  rate  of  increase  of  the  contents  of  the  tension- 
water  element:  the  instantaneous  moisture  input  and  the  evapo- 
transpiration  demand.  The  evapot ranspi ra t i on  rate  from  the 
upper-zone  tension-water  element  is  U2(x^/x^).  The  moisture 
input  rate,  up  to  the  point  when  tension-water  requirements 
are  met,  is  u-^.  At  the  time  when  x^  =  x° ,  two  possibilities 
arise : 

•  The  moisture  input  rate  is  larger  than 

or  equal  to  the  evapot ranspi ration  demand 

•  The  moisture  input  rate  is  smaller  than 
the  evapotranspi rat  ion  demand. 

In  the  first  case,  part  of  the  input  moisture  replaces  the 
amount  of  water  which  is  being  lost  through  the  evapot rans- 
pirafion  process.  The  tension-water  content  remains  at  its 
maximum  value  and  the  free-water  element  receives  the  excess 
moisture  input  (see  discussion  below).  In  the  second  ase , 
the  net  rate  of  increase  of  the  t ens i on -wa t e r  content  is  nega¬ 
tive  and  Xj  starts  to  diminish.  Summarizing, 
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U2 


|he<xl 


x°)  +  "  x°)  hf(u2  -  uj>] 

( A-  1  ) 


where  hg  and  are  given  by  Eqs.  2.2-23  and  2.2-24. 


The  moisture  input  rate  available  for  the  upper  zone 
free  water  element  is  (u^  -  u2>  h^(x^  '  x°)  *1e^u2  "  ui  )  •  At 
any  given  time  the  upper-zone  free-water  element  loses  water 
through  interflow  at  a  rate  of  dux2  and  supplies  the  lower 
zone  through  percolation  at  a  rate  p  given  by  Eq .  2.2-13.  If 
the  net  inflow  rate  to  the  upper-zone  free-water  element  is 
positive  for  a  prolonged  period  of  time,  x2  will  attain  its 
maximum  value  and  surface  runoff  will  start  to  occur.  The 
rate  at  which  surface  runoff  occurs  is  the  excess  of  the  mois¬ 
ture  input,  <Uj  -  u2)  hf(x  -  x°)  hg(u2  "  ) *  over  the  sum  of 

interflow  and  percolation  rates,  dux2  +  p.  When  this  quantity 
is  negative,  the  upper-zone  free-water  content  starts  to  diminish. 
Thus  , 


x2  =  ^U1  '  u2)  hf^xl  '  xl^  he(u2  ’  ul^  ‘  dux2  ‘  p* 

*  ^he^x2  “  x2  ^  +  hf<x2  ‘  x2  ^ 

x  hf|dux2  +  p  -  ( u  ^  -  u  2 )  h{.(Xj  ~  x°)  he(u9  -  Uj)| 

<A-2) 

Recall  the  constraint  on  the  ratio  of  free  to  tension 
water:  the  normalized  free-water  content  should  not  exceed 

the  normalized  tension-water  content,  i.e., 


.  o  .  o 

x  2  /  x  2  1  x  j  /  x  i 


(A- 3) 
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Therefore,  Inequality  A-5  is  automatically  satisfied. 

Consider,  next,  the  second  al ternat ive .  The  state 
equations  for  the  upper-zone  elements  become 


Xi<t)  =  Uj  -  u^x^/x^ 

*•>(  t>  =  -dx,  -  p 


( A-9 ) 

( A- 10 ) 


Therefore,  the  normalized  inflow  rates  are 


Xj(t)  u1  u2 

o  ~  o  ^  o 

X'i  Xj  Xj 


(A-n  > 


x2(t ) 


-  d  ill  - 
u  o 

X2 


( A- 12 ) 


In  this  case,  Inequality  A-5  is  equivalent  to 


■d 

u 


u,  u  0 

.  R_  <  _1  .  _2 

.o  —  .o  o 

*2  x  ^  x j 


(A- 13) 


which  can  be  rewritten  as 


u0  u , 

i  -—)*<—  +  £- 

U  0-0  o 

Xj  Xj  x2 


(A- 10) 


Since  the  expression  on  the  right  hand  side  of  last  inequality 
is  nonnegative,  the  inequality  will  certainly  hold  if 


d  >  u0/x, 
u  -  2  1 


(A- 15) 


Table  A-l  compares  the  values  of  the  interflow  param¬ 
eter  d^  with  the  maximum  evapot ranspi rat  ion  rate,  u2  ,  normal¬ 
ized  by  the  capacity  of  the  upper-zone  tension  water  element  , 
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TABLE  A- 1 

COMPARISON  OF  E VAPOTRAN S P I RAT 1 ON  AND  DRAINAGE  RATES 


BASIN 

o 

u2 

o 

X1 

0.0 

u2/xl 

du 

(mm/hr ) 

(mm) 

(1/hr) 

(1/hr) 

White  River 

5 . 8333E-2 

50 

1 . 1667E-3 

1 .A861E-2 

French  Broad 

5 . 0000E-2 

85 

5 . 882AE-A 

1 . A861E-2 

Bird  Creek 

2 . 7850E-2 

120 

2.3958E-4 

... 

1 .A861E-2 

Leaf  River 

5 . 6250E-2 

1 .79A9E-2 

5 . 73A7E-A 

1 .7313E-2 

Danville,  VA 

1 . 9250E- 1 

2A9 

7 . 7309E-A 

A.619AE-2 

Ariton,  Ala. 

1 .9167E-1 

75 

9.2976E-3 

Fulton,  Miss. 

2 . 2917E-1 

70 

3 . 2739E- 3 

3.3271E-2 

Culloden,  Ga . 

1 .9583E-1 

132 

1 . A836E-3 

1 . 79A9E-2 

3 . 3921E-3 

3 . 5  5  5  5  K  -  2 

1  .A583E-1 

10 

1 . A583E-2 

1 . A861E-2 

Eagle  River 

1 .8333E-1 

20 

9 . 1 665E-3 

1 . A861E-2 

S.  Yamhill  River 

1 . 1875E-1 

120 

9 . 8958E- A 

1 . 79A9E-2 

Clear  Boggy  Creek 

A.  1667E-2 

25 

1 . 6067E- 3 

3 . 6103E-3 

Illinois  River 

8 . 3333E-2 

28 

2.9702F.-3 

8.  1  K.GE-3 

Beaver  Creek 

5.A167E-2 

2 . 0002 E- 3 

1 . 3A57E-2 

Baron  Fork 

8 . 3333E-2 

21 

3 . 9682E-3 

1 . 1 107E-2 

Shoal  Creek 

6 . 2500E-2 

20 

3. 1250E-3 

1  .  5521  F.-2 

,  for  17  different  basins.  In  all  instances  Inequality  A- IS 
is  satisfied. 
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Finally,  for  the  last  possibility,  x^=x°  and  x2=x2‘ 

The  state  equations  become 

Xj(t)  =  ( u ^  -  u2)  hf(u2  *  Uj) 
x2 ( t )  =  ( (u^  -  u2)  he(u2  -  )  -  dux2  -  p) 

x  hf{dux2  +  p  -  ( u^  -  u  2 )  hg(u2  -  )  ] 

If  Uj  >  u2 ,  Eqs.  A-16  and  A-17  reduce  to 

x,(t)  =  0  ( A- 18 ) 

x2(t)  =  ( -  u2  *  dux2  *  P^  hf^dux2  +  P  "  ui  +  u2  ^ 

( A- 19 ) 

Thus,  x2(t)  £  0  and,  consequently.  Inequality  A-5  must  be  satis 
fied.  If  u^  ±  u2 ,  the  state-equations  ar.  identical  to  Eqs. 

A-9  and  A-10  and  the  same  analysis  used  for  the  case  0  <  ip  <  1 
shows  that  the  normalized  inflow  rates  into  the  upper-zone 
elements  satisfy  Inequality  A-5. 

Therefore,  for  all  basins  listed  in  Table  A-l,  the 
constraint  on  the  ratio  of  free  to  tension-water  for  the  upper- 
zone  is  ineffective. 


(A-16) 

(A-17) 
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